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INTRODUCTION 

Acoustic  or  sonic  fatigue,  the  deterioration  of  material  and 
structural  strength  from  high  frequency  wide  band  noise,  is  an 
important  factor  in  the  structural  design  and  safety  of  modern 
aircraft  and  aerospace  structures.  The  need  to  understand  how 
noise  is  generated  and  how  structures  respond  to  that  noise  is 
critical  to  the  design  of  future  aircraft  such  as  the  National 
Aerospace  Plane  (HASP)  and  various  short  take-off  and  landing 
aircraft  (STOL).  Advanced  engine  designs  using  propfans,  pres¬ 
ently  noisier  than  conventional  turbofans,  add  to  the  difficulty 
of  properly  designing  structures  without  the  support  of  an  ade¬ 
quate  experimental  database. 

The  need  for  predictive  techniques  that  are  accurate  and 
easy  to  use  is  of  paramount  importance  to  designers  of  present 
and  future  aircraft  structures.  Notwithstanding  the  uncertain¬ 
ties  associated  with  acoustic  load  generation,  stress  prediction, 
adverse  thermal  environments,  and  the  lack  of  fatigue  damage  pre¬ 
diction  models  for  both  metals  and  composite  materials,  it  is 
important  ^o  have  mathematical  models  that  can  represent  the 
random  nature  of  the  noise  environment  to  a  degree  that  permits 
rational  design  to  proceed.  Of  associated  concern  is  that  with 
increased  noise  levels,  the  structural  panels  can  behave  in  a 
nonlinear  fashion;  thus,  linear  mathematical  models  may  be  inap¬ 
propriate  for  use.  Indeed,  even  if  linear  models  were  conserva¬ 
tive  in  all  respects,  their  use  would  only  lead  to  designs  which 
were  inherently  too  heavy,  thereby  potentially  provoking  design 
modifications  that  are  more  costly  and  less  optimal. 

New  materials  such  as  polymer  and  metal  matrix  composites 
have  shown  themselves  to  be  potentially  useful  in  aircraft  struc¬ 
tures;  however,  the  limited  amount  of  data  related  to  their  per¬ 
formance  in  the  presence  of  an  acoustic  environment  is  disturbing 
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to  a  designer.  Thus,  analytical  techniques  for  predicting  life 
and  reliability  when  only  a  minimal  amount  of  material  property 
data  is  available  are  urgently  needed. 

The  surface  protection  systems  of  aerospace  and  aircraft 
structures  arc  usually  constructed  from  discretely  stiffened 
panels  or  stiffened  shells.  High  cycle  fatigue  failures  have 
occurred  in  these  structures  with  the  majority  of  fatigue  cracks 
appearing  in  the  near  vicinity  of  the  stiffening  element  or  the 
stiffener  Itself  (References  1-5).  Proper  dynamic  interaction 
between  the  panel  and  various  stiffening  elements  snould  be  taken 
into  account  when  calculating  the  response  of  the  panels  and  of 
the  stiffeners. 

Extensive  research  has  been  carried  out  for  linear  and  non¬ 
linear  analysis  of  a  single  bay  panel  (References  6-13).  For 
discretely  stiffened  panels,  most  analytical  work  is  based  on 
linear  theory.  However,  under  intensive  acoustic,  aerodynamic, 
and  thermal  loadings,  these  panels  vibrate  in  a  nonlinear  fashion 
and  a  nonlinear  analysis  is  needed  to  predict  deformations, 
stresses,  and  fatigue  life.  Different  methods  have  been  proposed 
to  study  random  vibrations  of  nonlinear  systems  (References 
14,15).  Among  the  most  widely  used  are  the  Fokker-Planck  equa¬ 
tion  solution  (References  16,17),  perturbation  method  (References 
18,19),  stochastic  linearization  (References  20-22),  and  the  time 
domain  Monte  Carlo  approach  (References  23-25).  Exact  solutions 
to  the  Fokker-Planck  equations  are  available  only  for  a  few 
simple  cases.  The  perturbation  method  is  usually  limited  to  one- 
or  two-degrees-of-f reedom  systems  and  is  valid  only  for  weakly 
nonlinear  cases.  The  stochastic  linearization  method,  although 
suitable  for  problems  with  strong  nonlinearities,  may  not  yield 
meaningful  results  for  complex  nonlinear  problems  that  might  be 
encountered  in  flight  structures.  The  time  domain  Monte  Carlo 
approach  can  be  used  efficiently  for  response  analysis  of  nonlin¬ 
ear  structures  subjected  to  random  pressure  fields*  In  this 
approach,  the  random  pressure  inputs  are  simulated  first  in  time 
domain  using  simulation  procedures  of  stationary  and  Gaussian 
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random  process  (Reference  23) »  then  the  resulting  nonlinear  equa¬ 
tions  of  notion  are  solved  in  the  tine  domain  by  numerical  tech¬ 
niques  (References  23-25). 

When  surface  protection  systems  of  a  flight  vehicle  are 
exposed  to  high  speed  aerodynamic  surface  flow  or  engine  exhaust 
hot  gases t  the  surface  temperatures  could  reach  3,000*F  (Refer¬ 
ence  26).  The  effects  of  those  high  thermal  gradients  are  degra¬ 
dation  of  strength!  stiffness,  and  fatigue  life  (References 
8,27).  In  addition,  structural-aerodynamic  instabilities  such  as 
buckling,  *'oil  canning,"  and  "snap  through"  could  be  induced  by 
the  action  of  thermal,  aerodynamic,  and  acoustic  loads  (Refer¬ 
ences  5,8,12,13).  These  effects  should  be  accounted  for  in  the 
nonlinear  response  of  structural  panels. 

The  objective  of  the  Phase  I  research  documented  herein  was 
to  develop  a  computational  procedure  for  predicting  the  life  and 
reliability  of  metal  and  composite  structural  panels  subjected  to 
complex  dynamic  loads  from  acoustic,  aerodynamic,  and  thermal 
environments.  Starting  with  existing  work  for  the  prediction  of 
sonic  fatigue  in  stiffened  pam Is  using  the  time  domain  method 
(Reference  28),  this  research  has  extended  the  model  to  include 
composite  materials.  Furthermore,  the  analytical  procedure  now 
includes  programs  to  compute  the  time  history  response  to  random 
noise  and  statistical  analysis  to  compute  probability  density  and 
peak  distribution  histograms  of  the  stresses,  uperossing  rates, 
and  expected  fatigue  damage. 
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II 


TECHNICAL  DISCUSSION 

The  basic  approach  used  in  the  present  research  can  be  cate¬ 
gorized  into  three  distinct  stages,  each  involving  its  own  set  of 
specialized  formulations.  The  first  stage  (Section  II  A.)  con¬ 
cerns  a  mathematical  description  and  phenomenological  representa¬ 
tion  of  the  acoustic  pressures  which  act  on  surrounding  structure 
(a  panel  in  this  work  is  used).  The  second  important  stage  (Sec¬ 
tions  II  B.  and  II  C.)  deals  with  the  kinematic  and  structural 
response  of  the  panel  to  the  acoustic  loading.  In  other  words, 
what  are  the  displacements  and  stresses?  The  third  stage  (Sec¬ 
tion  II  D. )  requires  an  estimation  of  panel  fatigue  life  based  on 
panel  structural  response  and  the  expected  acoustic  environment. 

The  relevant  equations  resulting  from  consideration  of  these 
three  stages  .are  contained  in  the  computer  code  discussed  in 
Section  II  E. 


A.  SIMULATION  OF  FJ^DOM  INPUT  PRESSURES  IN  SPACE>TIHE  DOMAIN 

The  random  pressure  acting  on  a  thermal  protection  system  of 
a  supersonic/hypersonic  aircraft  or  sub-orbital/orbital  vehicle 
arises  from  engine  exhaust  noise,  turbulent  surface  flow,  oscil¬ 
lating  shocks,  and  flow  separations.  In  addition,  there  might  be 
induced  structure-borne  dynamic  loads . due  to  engine  and  equipment 
vibrations  and  nonsteady  aerodynamic  loads  due  to  high  speed 
surface  flow  and  vibrations  of  the  thermal  protection  systems. 

During  a  normal  mission  consisting  of  take-off,  maneuvers, 
steady  flight,  and  landing,  the  thermal  protection  systems  will 
be  exposed  not  only  to  long  duration  stationary  pressure  but  also 
to  short-burst  intense  nonstationary  pressures.  Furthermore, 
oscillating  surface  shocks  and  flow  separations  could  create 
severe  localized  pressures  which  must  be  accounted  for  when  pre¬ 
dicting  the  fatigue  life  of  a  structural  component. 
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In  the  Phase  I  work  documented  herein •  the  analysis  of  the 
random  input  pressure  will  be  limited  to  stationary  and  Gaussian 
random  pressures  arising  from  engine  exhaust  noise  and  turbulent 
boundary  layer  flow.  Examination  of  the  nonstationary  and  non- 
Gaussian  characteristics  of  surface  pressures  arising  in  regions 
of  oscillating  shocks,  separated  flows,  and  rapid  changes  in 
thrust  requirements  will  be  addressed  during  the  Phase  II 
activities. 


1.  Simulation  of  Stationary-Homogeneous  Gaussian  Random 
Pressure 

Consider  a  random  pressure  pCx.y.O  acting  on  the  surface 
of  a  high  speed  flight  vehicle.  The  pressure  acting  normal  to 
the  surface  varies  randomly  in  time  and  space  along  the  surface 
coordinates  x  and  y.  The  pressure  p(5c,  y,0  is  characterized  by  a 
cross-spectral  density  function  5,(5, ti, a>)  where  and 

T^-yj-yj  are  the  spatial  separations  and  u>  is  frequency.  The 
cross-spectral  density  S,  is  obtained  utilizing  experimental 
data,  and  various  empirical  forms  are  available  for  jet  noise 
(References  5,29),  rocket  noise  (Reference  30),  and  turbulent 
boundary  layer  flow  (References  7,9).  The  simplest  form  of  the 
cross-spectral  density  is  the  truncated  Gaussian  white  noise 
pressure  uniformly  distributed  with  spatial  coordinates  x  and  y 
and 


11,03) 


if 

if 


0  2S03  ^  0)y 

03  <0  or  03>03„ 


(1) 


where  5,  is  a  given  constant  and  is  the  upper  cut-off 
frequency. 

The  spectral  density  5p(Ar,,^2,u>)  can  be  obtained  in  wave¬ 
number-frequency  domain  by  taking  the  Fourier  transformation  of 
as 
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SpCfci.fca.o})*  J  J  Sp(|,Ti,cx>)e'^*‘^e“^**’’d^dTi  .  (2 


‘  '  C2nyJ  J  - - 

^  '  -•  -to 

Then,  the  random  pressure  pCx.y.t)  can  be  simulated  by  the  series 
(Reference  23) 

w,  Wj  w,  , 

PCX .  y ,  0  -  V2 II  X  [-SpC/c  ,cOr)AA:iAfc2Aaj]“ 

i-l  y-i  r-l 

cos(A:ijX  + fcgyy +  a)j.t  +  <j)j^r)  (3 

where  are  realized  values  of  independent  random  phase  angles 

uniformly  distributed  between  0  and  2n.  The  values  of  the 
spectra  are  selected  at 


k  }j  ^  fcjj"^iA/Cj 

1,2,..*,  A/ 

^2/  ""  ^21^  J  ^^2 

7-  l,2,...,iV2 

00^ «  03,  +  rA<JU 

r- 

where  the  wave  numbei  and  frequency  intervals  are 

A/b2'“  C^2u””^2i)/  ^  2 

Aoo  “  (a\~ 

in  which  the  sutacrlpts  u  and  I  indicate  the  lower  and  the  upper 
cut-off  values  of  the  wave  number  and  frequency,  respectively. 


For  a  generation  of  random  sample  functions  from  Equation 
3  with  the  spectral  density  function  close  to  the  one  specified 
and  a  large  number  of  spatial  (Af.y)  and  time  (t)  points,  Ni,  Nz, 
and  Nz  must  be  large,  and,  as  a  consequence,  a  significant  amount 
of  computer  time  is  required  to  achieve  a  simulation. 

To  reduce  the  computation  costs  and  improve  the  effi¬ 
ciency  of  simulation,  the  Fast  Fourier  Transform  (FFT)  technique 
can  be  utilized  (References  31,32).  Rewriting  Equation  3  in  the 
form 


p(x,y  ,0  -  Rq 


i<,-l  i<2*l 


,  I 

L  i-o 


I 

/-o 


Z 

r-O 


(6) 


and  evaluating  pix.y.t'i  at 


x^^mAXt 

m*0, 1  1 

y  =  nAy , 

0, 1 . Mg-  1 

t  -  qAt, 

Q"0,1  ...,  Ml  3  “  1 

(7) 

where  "Re"  indicates  the  real  part  of  Equation  6  and 

1 

^  Ijr  ^  li  *  ^  21  •  ^ 

In  Equations  6  and  7, 

M2^  2."*^  "‘V2^  2^  ^2’‘  2"“ 

M3-2'”"-V3.V3>A^3-2'‘=’  (9) 
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and 


Ax"  211/MiAA:,  •2-si/yj 

Ay  *  2n/ MzAkz  -  2ii/ Vg.^au 

At  *  2jl/M3Aa)  =»  2E/V3a)u 


with  V2“2"‘*’"*,  V3“2™*"*,  and  mi>ni,mi>n2,  nia  >  713  all 

being  positive  integers.  The  use  of  the  FFT  technique  can  be 
applied  to  Equation  6  directly  resulting  in  a  drastic  reduction 
of  computer  time. 

A  second  difficulty  arising  from  the  simulation  of  a 
three  dimensional  random  process t  as  shown  in  Equation  6,  is  the 
creation  of  a  large  number  of  complex  numbers  that  rjust  be  stored 
in  the  computer  for  a  standard  FFT  application.  For  example,  if 
A^i  -  N2  -  JV3  -  128,  v,-Va-V3-4,  M  j  -  M2  -  Ma  -  512,  a  three  dimensional 
field  of  512  by  512  by  512  complex  numbers  must  be  stored  in 
order  to  perform  a  three  dimensional  FFT  procedure.  This  storage 
requirement  is  significant  even  for  large  mainframe  computers. 
Thus,  simulations  of  random  processes  have  been  primarily  applied 
to  one  dimensional  and,  in  some  cases,  two  dimensional  uses.  If 
a  random  pressure  acting  on  a  structural  surface  can  be  assumed 
to  be  uniformly  distributed  in  space.  Equation  6  reduces  to  a  one 
dimensional  simulation 


M-l 


r-O 


r 


e  e 


where 


yi,-[2Sp((A),)Aa)]' 


(11) 


(12) 


8 


and  S„(u>)  is  the  power  spectral  density  of  random  surface  pres¬ 
sure.  In  order  that  random  pressures  can  be  simulated  either 
from  Equation  6  or  Equation  11 »  the  input  spectral  densities  need 
to  be  prescribed. 


2.  Turbulent  Boundary  Layer  Flow 

Convective  turbulent  flow  produces  random  pressure  fluc¬ 
tuations  that  act  on  the  surface  protection  system  of  all  flight 
vehicles.  A  considerable  amount  of  work,  both  theoretical  and 
experimental,  has  been  carried  out  on  this  subject  with  regard  to 
panel  response,  panel  flutter,  and  noise  transmission  (References 
7,9,10,24).  However,  for  high  speed  supersonic  and  hypersonic 
flow,  the  information  that  is  available  seems  to  be  very  limited, 
and  a  substantial  amount  of  work  will  be  needed  to  characterize 
random  pressures  for  high  speed  flows.  For  the  purpose  of  this 
work,  the  semi-empirical  forms  of  the  cross-spectral  density  cor¬ 
responding  to  separated  supersonic  flow  given  in  Reference  9  will 
be  considered.  For  a  homogeneous  turbulent  flow  convecting  in 
the  x-direction  over  a  structural  surface,  the  cross-spectral 
denf ity  can  be  expressed  as 


(13) 

in  which  /? Cl . h. tw)*  *5(00),  and  £/*  the  correlation  coefficient, 
surface  pressure  spectral  density,  and  convection  velocity, 
respectively.  For  separated  supersonic  flow,  the  empirical 
formulas  from  Referer  9  are 
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i?y(0,Tl  ,CJO)  «  G 


(15) 


-«2lnl 


S(oo)-(6q:/^«)G 


(-8.094-1. 239 t»-0.295a)“-0.090a)®-0.0i4a>*-0.00liS’) 


(16) 


where  oo -  ln(a>6/2n£/.),  6,  £/«,  and  are  boundary  layer  thickness, 

free  stt  am  velocity,  and  free  stream  dynamic  pressure,  respec¬ 
tively.  The  attenuation  coefficients  ot,  and  aa  indicate  the 
degree  of  spatial  correlation  of  random  pressure  p(x,y,0* 

Similar  forms  of  cross-spectral  density  are  available  for  sub¬ 
sonic  and  attached  supersonic  flow  (References  7,9). 

For  a  linear  response  analysis  of  surface  panels  using 
the  power  spectral  density  approach,  Equation  13  can  be  used 
directly  as  an  input  parameter.  This  pi'ocedure  is  described  in 
Section  IT  C.  1.  However,  for  the  time-domain  nonlinear  response 
study  simulated  space-time  histories  are  needed  as  presented  in 
Equations  6  and  11.  The  required  spectral  density  can  be 
obtained  by  substituting  Elquation  13  into  Equation  2  and  perft  rni- 
ing  Fourier  transformation 


Sp(A:i  ,A:2,a)) 


S(a))ai  ‘Ua 

3t^[a?  +  (oj/(7c->-  A:i)^][Qt2''‘ 


(17) 


For  low  supersonic  flow  «,t  Mach  number  =  2,  0-  0.91  in.,  a,- 
1.22,  ttj-  0.26,  and  i7j“0.75£/.  (Reference  9). 


3.  Jet  Engine  Exhaust  Noise 

One  of  the  primary  causes  of  fatigue  in  many  flight 
structures  is  a  result  of  acoustic  loads  generated  by  near-field 
Jet  exhaust  noise.  Various  empirical  forms  similar  to  those 
given  in  Equations  13-16  are  available  to  characterize  the  random 
pressure  due  to  jet  exhaust  noise.  Kcwever,  for  supersonic- 
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hypersonic  aircraft  such  as  the  NASPt  detailed  statistical  infor¬ 
mation  on  the  localized  pressures  does  not  seem  to  be  available 
at  the  present  time.  Preliminary  estimates  indicate  that  the 
local  noise  levels  will  be  very  large »  exceeding  180  decibels 
(Reference  26). 


4.  Uniform  Distribution  of  Random  Pressure 

Useful  approximations  can  be  obtained  by  assuming  the 
input  pressure  to  be  uniformly  distributed  over  the  surface  of  a 
structural  component.  Theni  the  cross-spectral  density  can  be 
approximated  using  band-limited  Gaussian  white  noise  conditions, 
such  as  the  one  given  in  Equation  1.  The  expression  for  So  can 
be  written  as 

(18) 

where  Po  is  the  reference  pressure,  po-2.9xiO'*  psi  (0.00002A///n“), 
and  SPL  is  the  sound  pressure  level  expressed  in  decibels. 


B.  STRUCTURAL  MODELING 

The  governing  equations  of  motion  for  an  orthotropic 
composite  are  given  by  (Reference  33) 


d'^U)  d^W  d^U) 

d^F  d^F  d^w  ^  ^ 

dy"  dx^  dx^  dxdy  dxdy 


£l 

ay 


w 

“a ' 
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•Ml^P^dx.y 


(19) 
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li*F  ,  ^  »'F  li*F 


+  V^iV^  + 


Jjl^Y-t 

\dxdyj  d 


^w  d^w 
dx^  dy^ 


(20) 


where  Du,  D12,  D22)  are  bending  stiffnesses;  an*  0:12 »  <^s2i 

and  0*6  are  membrane  compliances;  Ni  and  Ny  are  the  inplane  loads 
applied  at  the  boundaries;  F  is  the  Airy  stress  function;  N^,  Ny, 
Ml,  and  My  are  the  inplane  and  bending  thermal  load  terms, 
respectively;  P'^  is  the  random  input  pressure;  nip  is  the  mass 
per  unit  area;  h  is  the  plate  thickness;  and  c  is  the  damping 
coefficient . 

The  structural  mechanics  embodied  in  Equations  19  and  20  are 
a  result  of  using  von  Karman*s  equations  for  large  deflections 
(References  33-35)  of  elastic  isotropic  plates  as  modified  by 
Ambartsumyan  for  orthotropic  materials  (Reference  33). 

The  terms  Nl,  Ny,  Ml,  and  Mj  are  computed  from 
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A/2 

-A/2 


n: 


M: 


A/2 

/  TCx,y,z)dZ 

cx  22  A  J 

-A/2 

A/2 

~  TCx,y,z-)ZdZ 


A/2 

f  T(x,y,z)ZdZ 


(21) 


(22) 


(23) 


(24^ 
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in  which  an  and  a22  are  coefficients  of  thermal  expansion  and  T 

is  the  temperature  distribution  in  the  panel. 

The  membrane  inplane  forces  are  given  by 


N. 


N, 


N 


xy 


b^F 

'  by^ 

(25) 

(26) 

bx^ 

b^F 

(27) 

bxby 

such  that  inplane  equilibrium  requirements  are  identically  satis¬ 
fied.  Equation  19  expresses  th<^  dynamic  equilibrium  in  the 
direction  normal  to  the  panel,  while  Equation  20  is  the  compati¬ 
bility  condition  for  inplane  strains. 

The  panel  is  assumed  to  be  simply  supported  on  all  four 
edges.  Exact  boundary  conditions  for  the  Airy  stress  function  F 
are  very  complicated  and,  for  the  present  study,  the  inplane 
boundary  conditions  are  satisfied  on  the  average  (References 
10,35,36).  Thus, 


b  a 

// 


du 

bx 


dxdy -  0 


(28) 
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dydx -  0 
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The  terms  u  and  v  in  Equations  28  and  29  are  the  inplane  dis¬ 
placements  which  are  expressed  as 


da  1 

^  d^F  d^F  ^ 

\  \(dwy 

dx  ' 

<  d^y  d^x  j 

1 

roi 

(34) 

dv  1 

^  d^F  d^F 

1 

ay 

0  X  d  y  j 

(35) 

Equations  28  and  29  imply  no  inplane  stretching  of  the  panel 
edges  in  an  average  sense  and  they  correspond  to  the  inplane 
boundary  conditions  for  immovable  edges.  Equations  30  and  31 
specify  that  the  average  inplane  shear  forces  are  zero  at  the 
boundary. 


C.  PREDICTION  OF  RESPONSE  IN  SURFACE  PANELS 

In  order  to  assess  the  fatigue  life  and  estimate  the  reli¬ 
ability  of  thermal  surface  protection  systems •  dynamic  response 
in  the  form  of  deflection  and  stress  is  needed.  For  linear 
systems,  response  calculations  can  be  obtained  either  in  the  time 
or  frequency  domains.  The  power  spectral  density  (PSD)  method  is 
commonly  used  to  obtain  solutions  for  a  frequency  domain 
approach.  For  a  time  domain  analysis,  the  simulation  of  random 
input  pressure,  as  described  in  Section  II  A.,  and  numerical 
integration  procedures  must  be  utilized.  When  the  response  is 
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nonlinear  I  a  time  domain  Monte  Carlo  type  method  can  be  developed 
to  obtain  deflection  and  stress  response  solutions  (References 
11,24,25,28.36>42) . 

In  the  present  study •  the  time  domain  analysis  is  verified 
for  a  linear  case  by  a  direct  comparison  of  the  response  predic¬ 
tions  to  those  obtained  by  the  PSD  procedure.  Then,  a  detailed 
study  of  nonlinear  response  using  a  Monte  Carlo  time  domain 
approach  is  developed. 

In  addition  to  the  very  high  noiy »  levels  that  will  be  acting 
on  the  surface  of  the  NASP-type  veh^oiesi  the  surface  tempera¬ 
tures  are  expected  tc  exceed  3,Q00*P  (Reference  26).  Such  high 
temperatures  will  induce  large  thermal  stresses  and  instabilities 
(buckling)  of  the  surface  panels.  However,  the  thermal  effects 
and  inplane  loads  are  not  considered  in  Phase  I  work. 


1  .t  The  Power  Spectral  Density  Method 

Consider  a  rectangular  panel,  shown  in  Figure  1,  exposed 
to  a  random  pressure  p(x.y,t). 

For  a  homogeneous  panel,  the  governing  equation  of  motion 
for  small  deformations  can  be  written  as 

DpV*w+ fiw  + m.pW’^  p(x  ,y  ,t) 


where 


Di,-£j,/iV12(1-v,2V2,) 

^12  "  ^21  ^  U 

D22-£'2,/iV12(1-V^2V2,) 
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in  which  £  n,  £22*  ‘v*  1^*  And  m,  are  moduli  of  elaaticity« 

Poisson’s  ratio,  damping  coefficient,  and  mass  per  unit  area, 
respectively.  The  solution  for  panel  deflection  can  be  expressed 
as  a  superposition  of  <.>rthogonal  modes  as 


-U.y.O-II  ^  mnCO-^  mnC^  *  y)  (37) 

in»  In- I 

where  q^a  the  generalized  coordinates  and  X„ut  Are  the  modes. 
Taking  the  Fourier  transformation  of  Equations  36  and  37  and  uti¬ 
lizing  orthogonality,  it  n  be  shown  that 


Q  mn  ^  mn  •*  mn 

where  the  frequency  response  function 

and  the  generalized  random  forces  are 


(38) 


(39) 


P 

in  which  a  bar  denotes  a  transformed  quantity.  The  modal  damping 
coefficients  And  the  natural  frequencies  of  panel  vibrations 
can  be  determined  from 


mn 


-7- f  f  p{x,y,^)X^^{x,y)dxdy 
m,oJ 0  Jo 
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where  is  a  modal  damping  coefficient  (percent  of  critical 

damping)  and  y  is  a  parameter  based  on  experimental  data#  For  a 
stationary  and  Gaussian  random  pressure  input  pC^.y.Oi  the 
deflection  response  spectral  density  can  be  determined  from 
(Reference  43) 


S.Cx.y.oo)-  t  t 


m-l r-1 n-l  i-1 


(43) 


where 


j  rm  ra  rh  rt 

^  mnri  ""  2 

in 0*0  JO  JO  Jo 


(44) 


where  an  asterisk  indicates  a  conjugate  quantity. 

The  orthonormal  modes  for  a  simply  supported  panel  are 


2 

^  mn  C  ^  ♦  y)* ^^sin(mjix/a)sin(reiiy/b)  . 

For  a  clamped-clamped  panel,  a  rough  approximation  of  the  modes 
can  be  obtained  by  using  clamped'-clamped  beam  modes 

X^(.x,y'j-X^\x)X„\y) 

(46) 

where 
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The  values  of  the  constants  Am$  Ym«  can  be  obtained  from 

Table  1.  Theni  for  a  simply  supported  panel  and  uniformly  dis¬ 
tributed  random  pressure  where  Sp(|,Ti,tt>)  -  S((a>),  the  cross- 
spectral  densities  of  generalized  random  forces  are 


-  (- 1 )"][  1  -  1 )"][  1  -  (- 1 )'][  1  -  (- 1 )']  • 

mpTLmnrl 


(48) 


For  a  clamped-clamped  panel 


Stresses  in  a  thin  panel  undergoing  linear  deformation 
can  be  calculated  from 
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TABLE  1 

CONSTANTS  FOR  CLAMPED-PLATE  MOOES 


m 

Vm 

<»„ 

1 

0.7133 

4.730040 

0.132857 

0.982 

2 

0.7068 

7.853202 

-0.0278749 

1.000 

3 

0.7071 

10.99561^ 

-0.00579227 

1.000 

4 

0.7071 

14.137164 

0.0012041 

1.000 

5 

0.7071 

17.278768 

0.002603 

1.000 

>5 

0.7071 

Yi-*-  (/n-5)it 

vm 

Sin  — 

1.000 

sinh^ 

12Z3 


(50) 


12D. 

Oy  —3 — z{d^w/ dy*  ^^d*w/ dx*) 

f  C 


(51) 


\ZD  ^  2 

- TS — w/dxdy 


(62) 


where  0,  and  Oy  are  the  normal  atreaaea  and  x.y  la  the  shear 

stress;  z  is  the  distance  from  the  mid-plane  of  the  panel. 

Taking  the  Fourier  transformation  of  Equations  60-52 1  using  Equa¬ 
tion  37 1  it  can  be  shown  that  the  spectral  densities  of  the 
stress  components  at  the  surface  of  a  simply  supported  panel  are 
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Sa^(x,y,a>) 
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L  “  J  nv-In-l 

I  Il"m»l*C'nn/a)*(r»n/b)»[a*x^/dxay]*-5‘^  . 

.  •*  J  fn«ln>I 


(64) 


(55) 


Equations  53-55  were  obtained  under  a  condition  that  the  cross- 
modal  terms  can  be  neglected  and  the  cross-spectral  densitieu  of 
the  generalized  random  forces  are  determined  from  Equation  44  by 
setting  m-r  and  n“l. 

The  root-mean-square  (rms)  values  of  displacement  and 
stresses  can  be  calculated  from 


rms 


/ 
L  0 


S  f-y  iUo^dw 


(56) 


where  S*  is  the  response  spectral  density  given  in  Equation  43 
for  displacements  and  Equations  53-55  for  stresses. 


2*  The  Time  Domain  Method 

To  solve  Equations  19  and  20,  panel  deflections  are 
expanded  in  terms  of  panel  modes: 


m-i n-l 

in  which  A^n  the  modal  amplitudes  and  <!>„„  are  the  natural 
modes  corresponding  to  a  linear  panel.  For  a  simply  supported 
panel,  may  be  written  as 


TT\  ^1,2, 3  , . . . 

n"  1 ,2,3,... 


(58) 
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where  Gin(mnx/a)  and  -  sin(nny/6)  . 

Substituting  Equation  67  into  Equation  20  yields: 
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(59) 


The  solution  for  F  consists  of  homogeneous  F^  ,  and  particular 
solution!  F p  .  The  particular  solution  for  F  can  readily  be 
obtained  from  Equation  59  by  following  the  procedures  given  in 
References  10  and  35: 
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where 
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For  the  case  of  inplane  boundary  conditions  that  are 
satisfied  on  the  average f  the  homogeneous  solution  can  be  written 
as 


^  (62) 

where  C Cz,  and  C3  are  constants  of  integration  that  are  to  be 

determixied  from  the  boundary  conditions  specified  in  Equations 
28-31.  Using  Equations  28-31  and  34  and  35,  together  with  the 
expresiiion  for  Airy  stress  function  F  and  the  solution  for  panel 
w.  the  consts.r*ts  of  integration  are 


(63) 


Cl  022^-0,2^ 

o  m  V  \  b* 


(64) 

(65) 

Having  completely  determined  Ft  Equation  19  is  solved  using  the 
Galerkin  method  by  computing  the  integral  average  of  Equation  19 
weighted  by  each  term  of  Equation  57.  The  natural  frequencies  of 
the  panel  obtained  by  the  linear  theory  are  imposed  into  Equation 
19.  Since  are  the  natural  modes  of  the  panel i  then 


C3  =  0  . 


dx 


- +^22- - 


dy‘ 


^p^mn^mn  (66) 


in  which  00^  are  the  natural  frequencies  of  the  undamped  linear 

system.  Substituting  and  w  as  given  in  Equation  57  into 
Equation  19 1  and  utilizing  the  Galerkin  type  procedure!  the  fol¬ 
lowing  system  of  nonlinear  differential  equations  are  cbtained  ' 
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where  a>,/  are  the  natural  frequencies  of  a  rectangular  panel. 
The  generalized  mass  and  the  generalized  random  forces  are 


a  b 
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The  nonlinear  stiffness  coefficients  Zu,d/nuri  are  given  by 

iQdfmkrL  "*  C  ^  0  [  iqd/mkr  i  C  ^  +  iqdfmkrl^.^  “  T  ,  fc  “  i)] 

where 
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^  iqdfnUtrlC^  ♦  " 
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^  UjdfnUert 


C-//-0 

otherwise 
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-  1 ,  7  -  -A:  0 

0 ,  otherwise 


/2,  i-Zc-O 
vCy*^)-\  1 »  7  “A:  9*0 

VO,  otherwise 


Before  step-’by-step  numerical  integration  can  be  ioiplem- 
anted  for  the  coupled  system  of  nonlinear  differential  equations 
given  in  Equation  67,  the  time  histories  of  the  genex’al iaed 
random  forces  Qt/(t)  are  needed.  This  is  achieved  by  first  simu¬ 
lating  the  multidimensional  random  pressure  p’'ix,y,t')  in  space¬ 
time  domain  and  then  evaluating  Equation  69  numerically  for  each 
value  of  indices  l.J  .  Following  the  procedures  given  in  Refer¬ 
ences  23  and  31,  the  stationary  random  pressure  p'‘(.x,y,t)  can  be 
simulated  as  given  in  Equation  3. 

An  alternate  procedure  to  generate  the  generalized  random  forces 
QtjiO  is  to  substitute  Equation  3  into  Equation  69  and  evaluate 
all  integrals  in  closed  form.  The  generalized  random  forces  are 
then  simulated  as  multi-variate  random  processes  (References 
31,32).  If  the  pressure  distribution  is  uniform  over  the  panel 
surface,  the  simulation  procedure  reduces  to  a  single  variate  and 
one  dimensional  random  process. 

The  displacement  and  stress  response  time  histories  are 
developed  for  one  realization  of  the  random  input  pressure 
p'’(x,y,0*  These  solutions  would  need  to  be  repeated  for  a  number 
of  different  realizations  and  then  the  response  statistics  calcu¬ 
lated  using  ensemble  averag  s  in  a  Monte  Carlo  sense.  However, 
by  assuming  the  input  pressure  to  be  an  ergodic  random  process, 
it  is  sufficient  to  obtain  a  solution  for  only  one  realization 
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and  then  use  temporal  averaging  to  calculate  the  required 
response  statistics.  Then*  the  mean  and  the  rms  values  of  the 
displacement  can  be  determined  from 


r 

j  wix,y,t)dt 
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(71 


,y 

\  0 

where  T  is  the  period  of  the  simulated  time  history  of  input 
pressure.  Similar  expressions  can  be  developed  for  panel 
stresses.  The  root-mean-square  values  calculated  from  Equations 
56  and  72  should  yield  the  same  results. 
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D.  SONIC  FATIQUE  OF  SURFACS  PANELS 

The  key  elements  in  predicting  the  fatigue  life  of  a  struc¬ 
tural  component  to  random  loads  are  detailed  stress  load  spectra 
at  a  critical  point  on  the  structure  and  reliable  cumulative 
damage  rules  for  random  stress  amplitudes.  For  a  multidimen¬ 
sional  stress  state*  the  most  damaging  stress  components  must  be 
known  as  well  as  the  choice  between  the  nominal  stress  and  the 
actual  load  stress  in  complex  geometries  and  connections.  The 
load  spectra  is  a  function  of  mission  requirements*  flight  condi¬ 
tions,  and  flight  duration.  The  information  on  threshold  cross¬ 
ings  and  peak  exceedances  is  needed  for  the  development  of  stress 
load  spectra.  In  addition  to  this  information*  stress  data  in 
the  form  of  S-N  diagrams  are  required.  Such  data  are  usually 
obtained  from  coupon  testing  under  constant  amplitude  loading  and 
are  approximated  by  the  following  equation 
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(76) 


where  S  is  a  fixed  stress  amplitude  for  constant  alternating 
loadi  N  is  the  number  of  stress  cycles  until  failure  at  the 
stress  level  5.  and  \  and  B  are  material  constants  depending  on 
the  material. 

Since  stress  response  of  surface  panels  is  random •  fatigue 
data  from  random  tests  should  give  the  required  information  for 
life  predictions.  In  this  case,  the  stress  response  S  is  repre¬ 
sented  as  the  root-mean-square  (rms)  value  and  the  number  of 
cycles  N  correspond  to  cycles  of  a  frequency  at  a  dominant 
response  peak.  For  a  linear  and  narrow  band  Gaussian  response, 
reasonable  predictions  of  fatigue  life  can  be  expected  using  this 
approach.  However,  under  severe  acoustic  and  thermal  environ¬ 
ment,  the  stress  response  is  nonlinear  and  non-Gaussian.  Fur¬ 
thermore,  for  most  practical  applications,  fatigue  data  under 
actual  random  inputs  are  rarely  available  for  full  scale  struc¬ 
tures  or  structural  components.  Most  of  the  fatigue  data 
digested  into  the  S-N  diagram  form  are  for  coupon  specimens  under 
constant  amplitude  loading.  This  information,  together  with  the 
distribution  of  stress  response  peaks,  can  be  utilized  to  con¬ 
struct  a  life  prediction  model  for  random  stresses. 

Consider  that  the  number  of  fatigue  cycles  is  equal  to  the 
number  of  positive  stress  peaks,  or  stress  reversaiu,  and  that 
each  damage  occurs  at  each  positive  stress  peak.  Then,  according 
to  the  Palmgren-Miner  linear  cumulative  damage  rule  (Reference 
44),  the  total  cumulative  damage  can  be  written  as 

D-  In(5,)/W(5,) 

i  (77) 

where  n  is  the  number  of  stress  peaks  (cycles)  experienced  by  the 
structure  at  stress  level  S,  ,  and  is  the  number  of  cycles  at 
which  failure  occurs.  Fatigue  failure  occurs  when  D  reaches  a 
value  of  unity.  Substituting  Equation  76  into  Equation  77  gives 
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(78) 


To  account  for  positive  peaks  that  occur  in  the  negative  stress 
region  and  add  to  the  cunulative  damage,  Equation  78  can  be 
modified  to 

.  (,8) 

Since  stress  response  is  a  random  quantity,  the  expected  damage 
in  the  time  interval  x  can  be  written  as  (Reference  43) 

f[gCc)]--^--  j  lSl'-p,CS)dS  (80) 

where  £*[^147]  is  the  expected  total  number  of  positive  stress  peaks 

and  P/(S)  is  the  probability  density  function  of  the  peak  magni¬ 
tude  of  the  random  stress  process.  The  expected  total  number  of 
peaks  can  be  estimated  from  (Reference  43) 

•  0 

FCMj]-- J  J  sP.^^(s.0,s)d9ds  (81) 

where  i*.,,  is  the  Joint  probability  density  of  S(0  (stress),  -$(0 
(stress  velocity),  and  S(0  (stress  acceleration). 

1.  Linear  Stress  Response 

The  response  of  a  linear  system  to  a  Gaussian  input  is 
also  Gaussian.  For  a  stationary  Gaussian  random  process,  the 
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probability  density  P,  and  the  joint  densities  P,t  ,  are 
known.  Theni  the  expected  niunber  of  total  peaks  and  the  density 
of  peaks  can  be  evaluated  in  closed  form  (Reference  43) 


(82) 


where  a,,  a^,  o«  are  standard  deviations  of  SCO*  <^CO*  <^(0 
parameter  a  is 


If  the  spectral  density «  S,((a>),  of  the  linear  stress  process  S(0 
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The  standard  deviations  a,,  Ct,  and  cr«  correspond  to  the  principal 
(naxinua)  stress  in  the  panel.  Two  liaiting  cases  of  peak  dis¬ 
tribution  can  be  obtained  for  a-1  and  a-0.  A  value  of  a  -  1 
corresponds  to  a  narrow  band  process  and  the  peak  distribution 
reduces  to  the  well  known  Rayleigh  distribution  while  for  very 
saall  a.  Equation  83  aay  be  approxiaated  by  a  Gaussian 
distribution. 

For  the  case  of  narrow  band  stress  response  a  -  1 .  the 
expected  nuaber  of  peaks  per  unit  tiae  reduces  to 


where  F[A^.(0)]  is  the  expected  rate  of  upcrossing  of  stress 

process  5(0  at  zero  threshold  level.  The  distribution  of  stress 
peaks  is  that  of  Rayleigh  distribution 


PM 
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The  expected  damage  can  be  obtained  in  closed  form  from  Equations 
80,  86,  and  87  as 


T(>/2a.)*'  r\+2\ 

f[D(-c)]-£[W.(0)] . •-  r(^— J  (88) 

where 

r(y)-2J  y>0  (89) 

0 

is  the  Gamma  function.  For  a  single  degree  of  freedom  response 
the  expected  upcrossing  rate  ^[Ar.CO)]  can  be  replaced  with  the 
natural  frequency  /q  (cycles/sec)  of  the  surface  panel.  For  the 
cases  where  many  modes  participate,  the  single  mode  approximation 
is  not  valid  and  expected  damage  should  be  calculated  using  Equa¬ 
tions  80,  81,  and  83.  In  this  case,  it  is  difficult  to  obtain  a 
closed  form  solution  and  a  numerical  procedure  is  used  to  evalu¬ 
ate  the  required  integrals.  The  expected  total  damage  in  the 

interval  from  x-Oto  x-T  (a  selected  time  period)  can  be 
r 

obtained  from  / £'[i?(x)]cix  .  For  a  stationary  random  response 

0 

process,  FCDCT)] -  TF[£)Cx)]  . 

As  shown  in  Reference  39,  the  expected  damage  ^[DCT)] 
does  not  change  by  much  when  plotted  versus  rtMTlT*  for  a  ranging 
from  0.25  to  1.  However,  for  a  wide  band  process  the  total 

number  of  stress  peaks,  EiMrl  ,  is  much  larger  than  for  a  narrow 
band  process  (a-l).  Thus,  fatigue  time  to  failure  at  £’[£)(T)]- 1 
will  be  shorter  for  a  wide  band  response.  or  example,  for  an 
extreme  case  when  a  is  very  small,  distribution  of  peaks  may  be 
approximated  by  a  Gaussian  probability  density 
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-00  <  S  <  00  . 


(90) 


V2iio, 


Then,  from  Equations  80  and  88 
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Since  i?[Mr]»>  fCiV.CO)],,  (Og)„>(o,)„,  where  the  subscripts  w  and  n 

indicate  wide  band  (many  modes)  and  narrow  band  (single  mode) 
stress  response,  the  expected  damage  for  a  wide  band  response 
will  be  larger  than  for  a  narrow  band  response.  However,  except 
for  these  extreme  cases  the  peak  distribution  of  a  stationary 
Gaussian  random  stress  process  is  neither  Rayleigh  nor  Gaussian. 


2*  Quasi -Nonlinear  Single  Degree  of  Freedom  Stress  Response 
Approximate  solutions  for  expected  fatigue  damage  can  be 
developed  by  assuming  the  nonlinear  panel  response  to  be  domi¬ 
nated  by  a  single  mode,  reducing  the  nonlinear  equations  of  panel 
motions  to  a  Duffing’ s  type  equation  and  using  a  linear  stress- 
strain  relationship.  Such  a  procedure  might  not  produce  meaning¬ 
ful  fatigue  estimates  of  a  realistic  surface  panel,  but  it  could 
give  preliminary  guidelines  on  the  effect  of  nonlinearities.  If 
the  inplane  motions  of  a  panel  are  constrained  at  the  edges,  the 
governing  equation  might  be  approximated  by  (Reference  45) 

D^*w-mpAc^{d^w/ d^w/ dy^)  +  CiW 

+  mpUj^  pCx,y,t)  (92) 


where 
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a  b 

^"2^/  f 


(93) 


and  c  is  the  wave  speed 


p(l-v=^) 


(94) 


where  p  is  the  material  density  and  v  is  the  Poisson’s  ratio. 
For  simple  support  boundary  conditions  and  single  mode 
approximation 


w^x  ,y  tt)  •  Q(Osin(rt:x:/a)sin(3ty/6)  . 


(95) 


Substituting  Equation  95  into  Equations  92  and  94,  and  utilizing 
orthogonality  gives 


g  +  11^ 11  Yq^)  “  -P(0 


(96) 


where 


4  r  C  Jix  ny 
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(97) 


Y-3/2/1'*  . 


(98) 


For  a  uniform  noise  pressure  distribution  over  the  panel  surface. 


P(t) 


rripii' 
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(99) 
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Stresses  in  the  panel  (linear  stress-strain  assumption)  can  be 
calculated  from  Equations  50-52  in  terms  of  the  generalized  coor¬ 
dinate  qf(f).  At  z^*h/2  and  the  middle  of  the  panel  ^x-a/2, 

y-6/2j,  x«y“0  and  ^he  principal  stress  is  the  larger  of  and 
Oy  .  Then,  the  principal  stress  S(i)  can  be  related  to  q(t)  as 


(100) 


where 
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(101) 


It  should  be  noted  that  Equation  100  is  a  crude  linear  approxima¬ 
tion  to  relate  stress  and  nonlinear  displacement.  Then,  Equation 
96  can  be  written  as 


S  +  +  1/(0 


(102) 


where  1/(0  is  assumed  to  be  a  Gaussian  white  noise  in  which 


and 


/Tlplt 


2L 


(103) 


The  approximate  solution  for  probablility  density  of  peak  magni¬ 
tude  was  obtained  in  Reference  43 
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PjiS) 


(104) 


2^oCO  n 

nk 


(S  +  gS^)© 


and 


where  is  the  spectral  density  level  of  the  Gaussian  white 

noise  aproximation.  Since  the  standard  deviation  of  linear 
stress  response  is 
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Equation  104  can  be  written  as 
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Substituting  Equation  107  into  Equation  80 1  the  expected  damage 
for  the  simplified  nonlinear  system  can  be  determined  from 


£[D(-c)]- 
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(108) 


The  expected  number  of  total  peaks  per  unit  time  for  a  narrow 
band  approximation  can  be  computed  from  (Reference  43) 


(109) 
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where  the  constant  C  can  be  evaluated  from  the  normalization 
condition 


m  o» 

j  j  P siis ,  s^dsds  "  I 


(110) 


where  is  the  joint  density  function  of  the  nonlinear  stress  s 
and  stress  velocity  i 


P,i(s,s) 


(111) 


Substituting  Equation  111  into  Equation  110  and  integrating 
(Reference  43) 


Vi 


(112) 


where  is  Bessel’s  function  of  the  second  kind  with  imaginary 

arguments.  To  the  first  order  in  €,  the  equation  can  be  approxi 
mated  as 

^ .  ... 

2jtciu„aJ  V  4  V  (11 

Thus  • 


2it(l  -ejoj) 


(114) 


39 


For  linear  response  e-0,  and  F[Ar.(0)]- a>n/2n  . 

To  evaluate  the  expected  damage  from  Equation  108,  numer¬ 
ical  integration  procedures  can  be  utilized.  The  Equations  108 
and  109  correspond  to  nonlinear  deflections  of  the  panel  with 
linear  stress-strain  relationships  assumed.  If  the  stress-strain 
relationship  is  nonlinear  such  as 

n-2  (115) 

where  are  proportionality  constants,  difficulties  would  arise 

in  obtaining  the  joint  density  function  and  the  peak  density 
function  P/(s).  For  n-2,  closed  form  solutions  were  obtained  for 
Pti,  EiMr]  and  P/(s)  in  Reference  43. 

The  procedure  presented  in  this  section  can  be  applied  to 
estimate  the  fatigue  life  of  simple  narrow  band  single-degree-of- 
freedom  systems.  The  linear  and  the  nonlinear  response  of 
surface  protection  systems  to  exhaust  noise  and  supersonic- 
hypersonic  turbulent  flow  will  be  composed  of  many  modes,  and 
simplified  models  presented  could  lead  to  erroneous  fatigue  life 
estimates.  However,  these  procedures  could  serve  as  useful 
guidelines  for  the  more  realistic  fatigue  life  estimates  of 
multimodal  nonlinear  systems. 


3.  Nonlinear  Stress  Response 

For  a  nonlinear  random  stress  response,  the  amplitude 
distribution  is  non-Gaussian.  Furthermore,  the  spectral  density 
of  stress  response  exhibits  a  wide  band  characteristic  indicating 
that  a  large  number  of  modes  could  be  contributing  to  the 
response  process.  An  improved  damage  model  for  nonlinear  stress 
can  be  constructed  by  computing  the  histograms  of  stress  peak 
distribution  directly  from  stress  response  time  histories  and 
then  using  Equation  80  to  predict  fatigue  damage.  The  total 
number  of  peaks  per  unit  time  £'[Mr]  also  can  be  evaluated  numer- 
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ically  from  the  response  time  history.  Such  a  procedure  has  been 
used  in  References  38-42  to  determine  fatigue  life  of  stiffened 
titanium  panels  at  room  and  elevated  temperatures. 

The  histograms  of  peak  distribution  and  the  total  number 
of  peaks  could  account  for  the  stress  overloads,  pressurization 
loads,  persistent  stress  reversals  due  to  snap-through,  oil 
canning,  etc.  However,  the  fatigue  relationship  given  in  Equa¬ 
tion  76  corresponds  to  either  a  constant  amplitude  stress  or  a 
root-mean-square  stress  of  a  stationary  process.  Furthermore, 
the  linear  damage  superposition  from  Equation  77  might  not  be 
valid  for  the  severe  loading  conditions  to  be  encountered  by  the 
surface  protection  systems.  For  metal  materials  and  low  cycle 
fatigue,  significant  improvements  have  been  made  in  predicting 
the  life  of  a  structure  by  utilizing  fracture  mechanics  and  time 
domain  stress  solutions  (Reference  46).  However,  for  high  cycle 
fatigue,  the  crack  propagation  stage  might  be  relatively  short  as 
compared  to  crack  initiation,  and  the  time  domain  crack  growth 
solution  might  not  be  meaningful  in  assessing  the  life  of  a 
structure.  Additionally,  for  composite  materials  no  reliable 
theory  seems  to  exist  for  predicting  crack  growth  as  a  function 
of  random  stress  response. 


E.  COMPUTER  PROGRAM  DESCRIPTION 

The  primary  objective  of  this  research  has  been  to  demon¬ 
strate  that  the  use  of  a  Monte  Carlo  simulation  of  a  random 
process  can  be  integrated  into  dynamical  equations  of  motion  to 
produce  a  time  domain  response  predictive  approach  to  understand¬ 
ing  the  fatigue  of  panels  exposed  to  acoustic  and  sonic  noise. 

The  theoretical  equations  of  the  time  domain  approach  presented 
in  Section  II  C.  1.  have  been  programmed  using  FORTRAN,  and  the 
resulting  computational  procedure  is  called  TDR  (Time  Domain 
Response).  In  this  section,  a  brief  explanation  of  TDR  will  be 
presented.  Appendix  A  of  this  report  contains  a  FORTRAN  listing. 
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TDR  is  written  to  handle  the  time  domain  response  analysis  of 
a  simply  supported  rectangular  panel  subjected  to  uniform  random 
pressure.  Once  TDR  has  been  successfully  executed,  an  output 
file  (FOR001.DAT)  is  created  which  contains  the  time  response 
history  of  displacements  and  stresses.  That  file  becomes  input 
to  another  program,  PDF,  which  calculates  the  probability  density 
functions  of  displacement  and  stresses  and  up-crossing  rates.  A 
schematic  of  these  programs  and  their  relationship  to  each  other 
is  provided  in  Figure  2.  The  files  created  by  PDF  can  be  read  by 
various  graphical  display  devices  to  produce  plots  of  statistical 
quantities  of  interest.  The  fatigue  life  of  a  panel  is  computed 
with  the  computer  program  DAMAGEl  which  uses  as  input  the  output 
file  FOR008.DAT  from  PDF. 
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RESULTS^ 

Earlier  sections  of  this  report  have  discussed  noise 
sources,  their  interaction  with  structure,  especially  light¬ 
weight  panels,  and  methods  for  determining  structural  response 
under  acoustic  and  sonic  noise  inputs.  The  primary  objective  of 
the  present  research  is  to  use  Monte  Carlo  simulation  techniques 
to  model  random  noise  in  such  a  manner  that  the  response  of  a 
panel  can  be  examined  in  the  time  domain.  The  validity  of  the 
time  domain  approach  is  verified  by  comparison  to  existing 
methods  and/or  experimental  data.  In  the  present  research  the 
lack  of  experimental  data  necessitates  the  former  manner  of 
verification. 

Since  existing  methods  of  sonic  fatigue  prediction  rely  upon 
linear  strain-displacement  relations,  the  theoretical  derivations 
in  the  time  domain  presented  earlier  were  reduced  to  equivalent 
linear  analysis.  The  results  from  the  linearized  time  domain 
method  is  then  directly  comparable  to  the  power  spectral  density 
approach.  The  results  of  that  comparison  are  provided  in  Section 
III  A.  which  is  concerned  with  isotropic  plates. 

A  significant  advantage  of  the  time  domain  method  presented 
herein  is  that  it  can  be  extended  to  the  regions  involving  both 
nonlinear  kinematic  and  nonlinear  material  behavior.  For  pur¬ 
poses  of  this  research,  nonlinear  kinematic  relations  were  used 
to  model  the-  strain-displacement  behavior  of  the  flat  panel.  The 
nonlinear  strain-displacement  relations  allow  the  modeling  of 
inplane  stretching  in  the  panel — a  phenomenon  which  is  extremely 
critical  to  the  prediction  of  fatigue  in  panels  exposed  to  high 
levels  of  acoustic  noise.  Demonstrations  of  the  importance  of 
this  modeling  assumption  are  provided  for  both  isotropic  and 
orthotropic  composite  panels  in  the  sections  to  follow. 
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The  extension  of  the  isotropic  derivations  to  the  composite 
idealization  is  important  and  underscores  the  versatility  and 
adaptability  of  the  time  domain  method.  While  not  addressed  in 
this  research t  the  role  of  transverse  shear  and  material  nonlin¬ 
earity  in  the  matrix  material  is  an  important  behavior  which  will 
be  examined  in  future  Phase  II  work. 

The  results  presented  herein  are  comprised  of  time  response 
histories  of  displacement  and  stress,  probability  density  and 
peak  distribution  histograms,  up-crossing  rates,  and  predictions 
of  sonic  fatigue  life  for  both  isotropic  and  oi  .hotropic  compos¬ 
ite  panels. 


A.  COMPARISONS  BETWEEN  THE  TIME  DOMAIN  AND  POWER  SPECTRAL 

DENSITY  METHODS  (ISOTROPIC  PLATES) 

As  an  example  of  the  veracity  of  the  time  domain  approach, 
the  nonlinear  response  and  fatigue  life  of  a  simply  supported 
panel  made  from  6A1-4V  titanium  material  is  predicted.  The  panel 
shown  in  Figure  1  is  assumed  to  be  exposed  to  a  uniformly  dis¬ 
tributed  stationary  Gaussian  random  pressure  for  which  the  trun¬ 
cated  spectral  density  is  given  in  Equation  18.  All  the  analyses 
presented  here  are  obtained  for  lower  and  upper  cut-off  frequen¬ 
cies  of  0  Hz  and  500  Hz,  respectively.  The  selected  frequency 
bandwidth  was  Au}-2n  rad/sec  (1  cycle/sec)  with  the  input  levels 
of  the  random  pressure  prescribed  in  decibels.  For  example,  a 
uniformly  distributed  white  noise  spectral  density  of  140  dB  over 
a  frequency  range  0-500  Hz  corresponds  to  an  overall  sound  pres¬ 
sure  level  of  167  dB.  If  the  upper  cut-off  frequency  is 
increased  to  1,000  Hz,  the  overall  sound  pressure  level  would 
increase  to  170  dB. 

The  random  input  pressure  p(0 '^as  simulated  from  Equation  11 
with  Ma  Af-4,096,  /Va-SlZ,  and  At -0.00025  second.  Thus,  the 
length  of  the  simulated  process  is  0.00025  second  x  4,096  =  1.024 
second.  The  fundamental  frequency  of  the  metallic  panel  selected 
for  this  study  is  about  100  Hz  and  the  fundamental  period  is  0.01 
second.  Thus,  the  simulated  process  covers  about  102  natural 
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periods  of  the  panel.  It  has  been  shown  in  previous  studies i 
that  for  a  stationary  response*  reasonable  statistical  properties 
can  be  obtained  fro*  a  tine  history  which  extends  for  about  100 
natural  periods  of  the  structure. 

The  numerical  results  were  obtained  for  a  panel  with  the  fol¬ 
lowing  geometric  and  material  properties:  a  =  20  inches*  b  «  8.2 
inches*  h  -  0.06  inches*  E  «  16.0  x  10*  psi*  v  *  0.34*  rnp'-pph  , 
pp  =  0.000414  Iby-sec^/in*.  The  modal  clamping  coefficients  were 
obtained  from  Equation  41  with  y  -  1  and  -  0.02.  The  numer¬ 
ical  calculations  were  obtained  for  three  modes  in  the  x  direc¬ 
tion  (/n- 1.2,3)  and  cne  mode  in  the  y  direction  (n- 1).  It  should 
be  noted  that  for  a  uniform  input  pressure  distribution  only  the 
odd  modes  contribute  to  panel  response.  The  modal  frequencies  of 
the  simply  supported  panel  are  CA}ii-620  rad/sec  (98.7  Hz)  and 
iUji  -  1 ,333  rad/sec  (212.2  Hz). 

When  the  panel  response  is  calculated  in  time  domain*  the 
spectral  densities  of  displacement  or  stress  process  Z(x*,y*, Oi 
where  x*  and  y*  are  the  selected  spatial  points  on  the  panel,  can 
be  obtained  utilizing  the  Fast  Fourier  Transiorm  (FFT)  technique 
(Reference  32).  The  finite  transform  of  Z(x*.y*,{„)  at  discrete 
frequencies  oo*  can  be  written  as  (Reference  32) 


in  which  T,  is  the  total  duration  of  the  response  random  process 

and  tn  -  nAt  vith  At  being  the  time  interval.  Then*  the  B'FT  numer¬ 
ical  estimate  of  the  response  spectral  density  is 
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The  numerical  results  presented  in  this  section  correspond  to  the 
center  of  the  panel,  i.e.,  x*  =  10  inches  and  y*  -  4.1  inches. 
Stresses  are  calculated  at  z-h/2  (top  surface  of  the  panel). 

The  displacement  response  time  histories  at  different  levels 
of  input  sound  pressure  are  shown  in  Figures  3  and  4.  For  an 
80  dB  (overall  s  107  dB)  pressure  input,  panel  response  is  linear 
and  the  largest  peaks  reach  about  0.01  inches.  At  120  dB  and 
higher  levels  of  input  pressure,  panel  response  is  nonlinear  with 
peaks  reaching  about  0.26  inch  for  150  dB  (overall  =  177  dB) 
input.  As  can  be  observed  from  these  results,  the  character  of 
the  random  response  process  changes  with  the  increasing  degree  of 
nonlinearity.  It  should  be  noted  that  the  mean  value  of  dis¬ 
placement  response  is  zero  for  the  linear  and  the  nonlinear 
cases. 

The  Oyy  stress  response  time  histories  are  presented  in 

Figures  5  and  6  for  several  different  levels  of  input  pressure. 
For  the  geometric  conditions  chosen  for  these  numerical  examples, 
the  shearing  stress  x^y^O  and  the  normal  o,.  stress  is  about  one 
half  the  value  of  the  Oy  stress.  Thus,  Oyy  is  the  principal 
stress. 

For  an  80  dB  sound  pressure  input,  panel  response  is  linear 
and  the  time  history  of  stress  response  is  similar  to  the  dis¬ 
placement  response.  However,  as  the  input  levels  increase  and 
the  panel  exhibits  nonlinear  vibrations,  the  stress  response 
changes  to  a  wide  band  process.  Furthermore,  the  mean  value  is 
not  zero  and  it  increases  with  the  increasing  sound  pressure 
input  level.  The  mean  values  of  stress  response  corresponding  to 
input  levels  80  dB,  120  dB,  140  dB,  and  150  dB  are,  respectively, 
0.04  psi,  431  psi,  5,418  psi,  anl  10,080  psi.  The  mean  stress  is 
caused  by  the  inplane  stretching  of  the  nonlinear  panel. 

The  nonlinear  transformation  between  stress  and  displacement 
(Equations  73,74,75)  contains  quadratic  terms  of  the  displacement 
component  w  ,  With  the  increasing  degree  of  nonlinearity,  these 
quadratic  terms  tend  to  dominate  the  stress-displacement  trans¬ 
formation  process.  These  effects  are  clearly  evident  in  Figures 
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5  and  6.  For  the  sound  pressure  input  of  ISO  dB  (177  dB 
overall))  the  stress  peaks  reach  67,000  psl»  It  should  be  noted 
that  even  these  high  stress  vslues  are  below  the  yield  strength 
(120)000  psi)  of  the  titaniuja  aaterial. 

The  root~iaean~aquare  (rms)  displacement  and  normal  stresses 
are  plotted  in  Figures  7  and  8  versus  the  rms  of  the  input  pres¬ 
sure.  The  linear  response  predictions  were  obtained  by  the  power 
spectral  density  (PSD)  and  the  time  domain  methods.  In  the  time 
domain  approach)  the  governing  equations  of  motion  and  the 
stress-displacement  relationships  were  reduced  to  those  of  a 
linear  problem.  However)  the  simulation  of  the  random  input 
pressure  and  the  solution  procedure  of  the  governing  linearized 
equations  were  identical  to  that  of  the  nonlinear  case.  As  can 
be  observed  from  these  results)  the  agreement  between  the  PSD 
method  and  the  time  domain  approach  is  very  good.  For  example) 
at  an  input  level  of  140  dB  (167  dB  overall)  the  rms  displace¬ 
ments  and  stresses  are:  Wrns  -  0.314,  inch)  0.306  inch;  o..  s 
14)30(/  psi,  13)900  psi;  s  26)760  psi)  26)020  pal  (PSD)  time 
domain) . 

It  should  be  noted  that  these  linear  response  predictions 
overestimate  the  actual  nonlinear  response  by  a  large  amount. 

This  can  be  seen  from  the  results  shown  in  Figures  7  and  8.  For 
the  input  levels  exceeding  about  110  dB  (137  dB  overall))  nonlin¬ 
ear  analysis  is  required  for  displacement  and  stress  response 
predictions.  This  input  limit  corresponds  only  for  the  simply 
supported  20  inch  by  8.2  inch  by  0.06  inch  titanium  panel  exposed 
to  uniform  random  noise  pressure.  For  stiffer  panels  this  limit¬ 
ing  input  pressure  value  would  increase  while  for  less  stiff 
panels  it  would  decrease. 

The  response  spectral  densities  for  deflection  and  normal 
stress  Oyy  are  shown  in  Figures  9-11.  At  low  input  levels  (80  dB 
and  100  dB)  where  response  is  linear)  distinct  peaks  can  be  seen 
at  the  modal  frequencies  of  98.7  Hz  and  212.2  Hz.  Furthermore) 
similar  characteristics  can  be  seen  between  the  spectral  densi¬ 
ties  of  displacement  and  stress.  As  the  input  pressure 
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increases t  the  distinct  peaks  that  are  characteristic  of  linear 
vibrations  tend  to  flatten  and  shift  towards  higher  frequencies. 
For  an  input  level  of  140  dBt  the  distinct  peaks  are  no  more 
evident  and  the  response  spectra  tends  to  exhibit  the  character¬ 
istics  of  a  wide  band  random  process.  In  addition •  the  shape  of 
the  displacement  and  stress  spectral  densities  is  different.  The 
spectral  densities  shown  in  these  figures  are  the  unsmoothed 
outputs  of  the  FFT  of  the  response  time  history  corresponding  to 
one  solution  realization. 

When  the  random  process  possesses  wide  band  characteristics, 
large  fluctuations  of  the  FFT  output  are  typical.  To  improve  the 
FFT  inverse  calculations,  a  larger  number  of  simulated  points  and 
smaller  time  steps  should  be  taken.  In  addition,  the  response 
spectral  densities  should  be  calculated  for  several  realizations 
of  the  random  input  pressure.  These  spectral  densities  are  then 
averaged  to  obtain  the  response  spectral  density. 

Displacement  or  stress  response  probability  density  histo¬ 
grams,  peak  distributions,  total  number  of  peaks  per  unit  time, 
and  thi'eshold  crossing  rates  can  be  obtained  from  the  response 
time  histories.  Figures  12  and  13  show  the  probability  density 
and  peak  distribution  histograms  for  the  nonlinear  displacement 
response.  For  comparison,  a  Qaussian  density  function  is  given 
with  each  probability  density  histogram  and  a  Rayleigh  distribu¬ 
tion  with  each  histogram  of  peak  distribution.  As  can  be 
observed  from  these  results,  the  nonlinear  response  is  no  longer 
Qaussian  and  the  peak  distribution  does  not  follow  the  Rayleigh 
distribution. 

Similar  results  are  presented  in  Figures  14  and  15  for  normal 
stress  component  .  For  the  nonlinear  stress  process,  large 
differences  can  be  seen  between  the  response  histograms  and  the 
theoretical  probability  and  peak  distributions.  These  large  dif¬ 
ferences  are  produced  by  the  nonlinear  relationship  between 
stress  and  displacement.  Thus,  the  various  approximate  theories 
which  predict  the  nonlinear  displacement  response  but  use  a 
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linear  atresa-diaplacement  relationship  to  obtain  streasea  do  not 
give  a  meaningful  procedure  to  obtain  nonlinear  stresses  for 
sonic  fatigue  analysis. 

From  the  results  presented  in  Figures  13  and  16,  it  can  be 
seen  that  for  a  highly  nonlinear  response  the  stress  process  con¬ 
tains  a  large  number  of  peaks  above  the  2o  (cr  =  standard  devi¬ 
ation)  range.  However,  the  histogram  of  displacement  peak 
distribution  does  not  show  any  peaks  above  the  2.2  cr  range.  Here 
the  peaks  seem  to  be  concentrated  at  about  the  1 . 5  o  value . 

Since  fatigue  life  is  very  sensitive  to  the  magnitude  of  stress 
ranges,  erroneous  fatigue  life  predictions  will  be  obtained  if 
the  form  of  the  stress  peak  distribution  is  assumed  to  be  the 
same  as  the  displacement  peak  distribution.  In  addition,  the 
nonlinear  stress  process  contains  a  mean  value  while  the  mean  is 
zero  for  a  nonlinear  displacement  response. 

The  threshold  up-crossing  rates  for  the  stress  process  are 

given  in  Figure  16  for  several  levels  of  input  pressure.  For  a 
linear  stress  response  at  80  dB  input,  the  threshold  up-crossing 
rate  closely  approximates  a  theoretical  Gaussian  prediction.  As 
the  input  pressure  increases  and  stress  response  becomes  more 
nonlinear,  the  up-crossing  rates  increase. 

The  expected  fatigue  damage  has  been  predicted  utilizing 
Equation  78.  The  expected  total  number  of  peaks  £"[ Mr]  were  esti¬ 
mated  directly  from  the  stress  response  time  histories.  The 
values  for  £'[Mr]  are  160,  228,  666,  and  748  peaks/second  for  120, 
130,  140,  and  160  dB  sound  pressure  Inputs,  respectively.  The 
integral  in  Equation  78  was  evaluated  numerically.  The  distribu¬ 
tion  of  stress  peaks  -P/Cs)  is  the  histogram  of  the  peak  distribu¬ 
tion  corresponding  to  the  principal  stress  (Figures  14  and 
15)  . 

To  Illustrate  the  fatigue  damage  estimation  procedure, 
typical  parameters  were  chosen  from  experimental  data  for  the 
titanium  material  under  room  temperature.  These  parameters  cor¬ 
respond  to  tests  at  constant  stress  amplitude  and  stress  ratio 
^--1.  The  stress  ratio  R  is  defined  as  in  which  Ob„„ 
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is  the  minimum  stress  and  is  the  maximum  stress.  A  stress 
ratio  corresponds  to  a  sero  mean  value.  The  experimental 

parameters  chosen  in  this  study  are  \-6.0and  5-1.518X10*.  When 
using  these  parameters •  the  stress  amplitude  is  in  units  of  ksi. 
To  account  for  stress  concentrations,  size  effects,  and  manufac¬ 
turing  imperfections,  a  fatigue  reduction  factor  /C/,  should  be 
introduced.  However,  for  the  simple  panel  chosen  in  this  study 
no  fatigue  reduction  factor  was  introduced. 

The  expected  fatigue  damage  is  plotted  in  Figure  17  versus 
the  product  of  FCMrJx  .  A  value  of  £f[Z?CT)]- 0.1  corresponds  to  10 
percent  damage  and  FCZ^Cx)]-  1  a  100  percent  damage  or  total 
failure  of  the  structure.  By  selecting  F[l>Cx)]-  1  and  the  calcu¬ 
lated  value  of  the  total  number  of  peaks  F[Afr]  «  the  time  to 
failure  can  be  obtained  from  Figure  17.  These  results  are  given 
in  Table  2 . 


TABLE  2 

STANDARD  DEVIATION  OF  cr^,  STRESS, 
EXPECTED  NUMBER  OF  PEAKS  AND  FATIGUE  LIFE 


Input  Sound 
Pressure,  dB 

Nusber  of 
Peaks/Second 

Standard 
Deviation,  psi 

Fatigue 

Life,  Hours 

Linear 

Nonlinear 

Linear 

Nonlinear 

Linear 

Nonlinear 

120  (147) 

124 

160 

2,676 

2,298 

2.39  X  10« 

938,888 

130  (167) 

124 

228 

8,462 

3,697 

2,393 

16,578 

140  (167) 

124 

566 

26,760 

9,260 

2.393 

18.7 

160  (177) 

124 

748 

84,622 

14,100 

— 

1.0 

(  )  s  Overall  levels 


When  the  stress  response  analysis  is  performed  using  a  linear 
plate  theory,  the  stresses  can  be  calculated  by  the  PSD  method 
and  the  expected  fatigue  damage  estimated  from  Equation  86.  The 
results  based  on  linear  theory  are  also  presented  in  Table  2. 
These  results  indicate  that  for  input  levels  above  120  dB  (147  dB 
overall)  the  linear  theory  over  estimates  stress  and  under  esti- 
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mates  fatigue  life  of  surface  panels.  For  sound  pressure  inputs 
of  120  dB  or  lower,  the  titanium  panel  would  not  fail  by  fatigue. 
Th  fatigue  life  estimates  presented  in  this  report  should  be 
viewed  as  preliminary  merely  to  illustrate  the  procedure  of  the 
time  domain  analysis. 


B.  PREDICTION  OF  RESPONSE  IN  ORTHOl’ROPIC  COMPOSITE  PANELS  USING 
lEE  TIME  DOMAIN  METHOD 

In  this  aection,  representative  examples  of  the  predictive 
capabilities  of  the  time  domain  approach  for  application  to 
orthotropic  panels  is  presented.  For  the  most  part,  the  examples 
provided  are  analogous  to  those  shown  for  the  isotropic  panel  of 
Section  III  A.  However,  examples  showing  the  effects  of  varia¬ 
tion  in  lay-up  angle  have  been  included  to  demonstrate  the 
design/analysis  capabilities  of  the  time  domain  approach. 

A  subtle  but  important  distinction  between  an  orthotropic  and 
laminated  composite  material  should  be  noted;  namely,  for  pur¬ 
poses  of  this  Phase  I  research,  the  orthotropic  derivation 
implicitly  assumes  that  the  entire  plate  is  composed  of  a  single 
material  possessing  orthotropic  properties,  as  contrasted  to  a 
laminated  panel  wherein  orthotropic  properties  generally  vary 
throughout  the  thickness.  This  assumption  presents  no  problem 
with  respect  to  the  calculation  of  displacement  response  time 
histories  and  statistics  since  the  constitutive  properties  of  the 
chosen  orthotropic  material  (ylu,  4i3i  Du,  Dix,  Dzz,  and 

have  been  obtained  from  the  equations  used  to  obtain  consti¬ 
tutive  properties  of  a  laminated  composite  structure. 

Derivation  of  stresses  within  an  orthotropic  or  laminated 
composite  panel  are  different,  however.  An  orthotropic  material 
has  a  membrane  stress  field  which  is  constant  through  the  thick¬ 
ness.  The  bending  stresses  vary  linearly  through  the  thickness 
with  the  stress  being  zero  at  the  location  of  the  neutral  axis. 
Thus,  for  an  orthotropic  material,  the  total  stress  is  simply  a 
function  of  the  through-thickness  coordinate  of  the  panel.  For  a 
laminated  composite  material,  the  stress  is  a  function  of  the 
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material  properties  of  the  particular  lamina.  Since)  in  general) 
lamina  lay-up  angles  vary  throughout  the  thickness)  a  laminated 
composite  requires  stress-displacement  relations  which  depend 
both  upon  the  through-thickness  location  of  the  lamina  and  the 
constitutive  properties  of  the  lamina. 

Figure  18  shows  an  example  of  the  stress  and  strain  distribu¬ 
tion  through  the  thickness  of  a  laminated  composite.  Note  that 
the  membrane  strain  field  is  constant  in  value  and  the  bending 
strains  vary  linearly  with  through-thickness  position.  The 
stresses  vary  according  to  the  lamina  constitutive  properties  and 
are  neither  constant  nor  linear.  The  stress  distribution  for  an 
orthotropic  material  would  be  similar  to  the  strain  distributions 
shown  in  Figure  18,  i.e.,  membrane  stress  is  constant  and  the 
bending  stress  is  linear  through  the  thickness. 

For  this  research,  the  simpler  orthotropic  stress- 
displacemen relntions  have  been  used  inasmuch  as  demonstration 
of  the  feasibility  of  the  time  domain  approach  was  the  primary 
objective.  During  Phase  II  the  relations  for  a  laminated  compos¬ 
ite  material  will  be  implemented. 

The  panel  shown  in  Figure  1  represented  the  basic  configura¬ 
tion  used  for  the  orthotropic  analysis  shown  herein.  As  before, 
the  lower  and  upper  cut-off  frequencies  are  0  and  500  Hz,  respec¬ 
tively.  The  duration  of  the  simulated  process  is  1.024  sec.  The 
first  laminated  composite  example  is  composed  of  the  lay-up 
[0/+ 45/ -  45/90],  for  a  total  of  eight  layers  and  an  overall  thick¬ 
ness  of  0.0416  inches.  Each  layer  is  made  frosB  A-S/3501  Graphi¬ 
te/Epoxy.  The  basic  lamina  constitutive  properties  are  obtained 
from  the  O-degree  lamina  elastic  properties,  i.e., 

El-  18.6F  +  06 
2.0E  +  06 
0.8F  +  06 
v«-  0.31 
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The  numerical  results  were  obtained  for  a  panel  with  overall 
dimensions  of  a  ~  20.0  inches*  b  =  8.2  inches,  »  0.0001302 
lb/-6ec‘/in^.  The  material  damping  factor  was  s  0.05.  The  cal¬ 
culated  modal  frequencies  were 


fu 

= 

133 

cycles/sec 

f  12 

= 

466 

cycles/sec 

/»> 

= 

1,021 

cycles/sec 

= 

200 

cycles/scc 

f  22 

533 

cycles/sec 

/  23 

1,088 

cycles/sec 

/  31 

= 

313 

cycles/aec 

/  32 

644 

cycles/sec 

/ss 

1,199 

cycles/sec 

As  an  example  of  the  importance  of  the  assumption  of  linear¬ 
ity  or  nonlinearity  with  respect  to  the  strain-displacement  rela¬ 
tions,  Figures  19  and  20  show  the  effect  of  the  two  assumptions. 
Figure  19,  the  linear  response,  shows  a  maximum  displacement  of 
almost  4  inches,  whereas  the  nonlinear  response  provides  only  0.3 
inches  approximately.  The  difference  is  dramatic  with  the  non¬ 
linear  response  obviously  the  more  realistic. 

In  the  x'esponse  histories  shown,  the  limitations  in  the 
graphical  display  device  allowed  only  240  discrete  points  of  dis¬ 
placement  and  time  to  be  plotted.  Thus,  a  separate  program  was 
written  to  take  the  time  history  produced  by  TOR,  containing 
4,096  discrete  points,  and  reduce  it  to  the  first  240  minimum  and 
maximum  points.  As  a  result,  the  time  histories  shown  appear 
slightly  different  than  those  prepared  for  the  isotropic  panel. 

In  comparing  Figures  19  and  20,  it  is  also  observed  that 
there  are  apparently  more  variations  between  high  and  low  values 
for  the  nonlinear  response  because  240  points  of  minimum/maximum 
appear  in  just  0.24  seconds  for  the  nonlinear  case  as  contrasted 
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to  about  0,63  seconds  for  the  linear  case.  The  Indication  that  a 
nonlinear  response  produces  more  changes  in  stress  direction 
could  have  important  implications  in  the  determination  of  panel 
fatigue  life  and  is  one  more  reason  why  the  nonlinear  assumption 
is  preferred  over  that  of  linear. 

The  time  histories  of  the  lateral  component  of  stress •  , 

for  both  the  linear  and  nonlinear  cases  are  shown  in  Figures  21 
and  22 t  respectively.  As  with  the  displacement  response,  the 
linear  prediction  provides  unrealistic  values,  whereas  the  non¬ 
linear  assumption  results  in  the  calculation  of  stress  values 
that  are  comparable  to  the  actual  expected  response  of  the  panel . 
Note  also  that  the  mean  value  for  the  nonlinear  assumption  is  not 
zero.  Clearly,  the  use  of  linear  strain-displacement  relations 
is  not  warranted  for  applications  involving  high  levels  of  acous¬ 
tic  loading. 

Displacement  and  stress  response  probability  density  histo¬ 
grams,  peak  distributions,  total  number  of  peaks  per  unit  time, 
and  threshold  crossing  rates  can  be  obtained  from  the  response 
time  histories.  Figures  23  and  24  show  the  probability  density 
histograms  for  two  different  sound  pressure  levels  (130  and  150 
dB)  using  nonlinear  strain-displacement  relations.  For  compari¬ 
son,  a  Gaussian  density  function  is  given  with  each  probability 
density  histogram.  The  nonlinear  response  is  no  longer  Gaussian. 

Figures  25  and  26  show  the  peak  distribution  histogram  for 
the  nonlinear  strain-displacement  relations  compared  to  a  Ray¬ 
leigh  distribution.  The  nonlinear  prediction  does  not  follow  the 
Rayleigh  distribution. 

Similar  histograms  of  probability  density  and  peak  distribu¬ 
tion  for  the  lateral  stress  component,  Oyy  ,  are  shown  in  Figures 
27-30.  As  for  the  isotropic  panel,  large  differences  between  the 
response  histograms  and  the  theoretical  probability  and  peak  dis¬ 
tributions  can  be  seen.  These  large  differences  are  produced  by 
the  nonlinear  relationship  between  stress  and  displacement  and 
are  additional  confirmation  that  the  assumption  of  linearity  in 


at rain-displacement  relations  is  not  appropriate  for  the  predic¬ 
tion  of  stress  response  in  panels  subject  to  high  levels  of 
acoustic  noise. 

Figure  31  shows  up-crossings  per  second  for  0yy  for  various 

levels  of  sound  pressure  level.  For  the  lower  sound  pressure 
level  (i.e.i  110  dB)  the  response  is  mostly  linear:  however,  for 
higher  levels  of  input  the  response  is  increasingly  nonlinear  and 
the  up-crossing  rate  increases  markedly. 

When  the  basic  lay-up  Just  used  is  varied  slightly,  it  is 
possible  to  alter  the  response  of  the  panel,  As  an  example. 
Figure  32  shows  the  Oyy  RMS  response  of  the  panel  for  150  dB  input 
when  the  lay-up  angle  of  the  interior  plies  are  varied  from  0 
through  90  degrees.  At  approximately  30  degrees  the  RMS  response 
is  at  a  minimum.  Figure  33  shows  how  the  up-crossing  rate  varies 
with  lay-up  angle — apparently,  the  up-crossing  rate  is  not 
strongly  dependent  on  the  lay-up  angle  for  this  particular  lami¬ 
nated  composite  example. 

As  a  final  example  of  the  usage  of  the  time  domain  approach, 
another  laminate  construction  is  examined.  Figures  34  and  35 
show  the  up-crossing  rate  for  stress  components  and  Oyy  for  a 
[+0/-03.a.  laminate.  The  effect  of  lay-up  angle  on  response  is  of 
importance  in  these  examples.  A  designer,  knowing  that  the 
luy-up  angle  would  have  an  effect  on  panel  response,  could  poten¬ 
tially  alter  the  design  in  such  a  manner  that  the  BMS  stresses 
and  up-crossing  rates  would  be  reduced  and  fatigue  life 
increased. 

Figure  36  shows  the  RMS  stresses  for  fiber  directed  and 
transverse  stresses  as  a  function  of  the  lay-up  angle.  The  fiber 
direction  stresses  are  relatively  low  in  this  example,  approxi¬ 
mately  15,000  psi,  when  compared  to  a  nominal  allowable  of 
180,000  psi.  However,  in  the  transverse  to  fiber  direction  the 
stress  level  peaks  at  about  6,300  psi,  which  is  very  dose  to  the 
static  stress  allowable  for  the  matrix  material.  In  an  actual 
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design  situation »  the  panel  would  need  sizing  to  ensure  that  the 
anticipated  stress  level  is  below  the  allowable  fatigue  stress 
level  for  the  matrix  material. 

Prediction  of  fatigue  life  in  composite  materials  is  not  pos¬ 
sible  at  this  time  because  of  limitations  in  the  theoretical 
understanding  of  fatigue  in  composites  and  the  lack  of  material 
property  data.  As  a  consequence •  predictions  of  sonic  fatigue 
life  for  a  composite  panel  have  not  been  produced.  Howeveri  the 
statistical  approach  documented  for  metal  panels  and  the  tech¬ 
niques  used  to  produce  predictions  based  on  time  domain  response 
analysis  are  valid  and  can  be  utilized  when  polymer-based  compos¬ 
ite  technology  advances.  Phase  II  work  will  develop  the  time 
domain  response  methods  for  laminated  composite  materials  to  the 
maximum  extent  possible.  Also,  Phase  II  work  will  concentrate  on 
the  appropriate  theoretical  relationship  to  be  used  to  predict 
composite  fatigue.  The  ability  to  predict  acoustically  generated 
stresses  combined  with  an  approach  for  predicting  fatigue  of  com¬ 
posite  materials  will  allow  an  analyst  the  ability  to  use  the 
time  domain  approach  to  predict  the  sonic  fatigue  life  of  struc¬ 
tural  panels i 
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IV 


CONCLUSIONS 

A  simple  rectangular  panel  was  selected  to  demonstrate  the 
applicability  of  time  domain  analysis  to  predict  nonlinear 
response  and  fatigue  life  of  metal  and  composite  panels.  It  is 
shown  that  the  linear  theory  overestimates  deflection  and  stress 
response  by  a  large  amount,  resulting  in  a  predicted  shorter 
fatigue  life.  If  one  were  to  use  linear  theory  as  a  design  tool, 
then  properly  designed  panels  would  be  relatively  stiffer  and 
heavier  than  required.  A  non-optimal,  weight-inefficient  struc¬ 
ture  would  result. 

The  nonlinear  response,  as  predicted  by  the  Monte  Carlo  time 
domain  approach,  is  non-Gaussian  and  peaks  do  not  follow  a  Ray¬ 
leigh  distribution.  With  a  nonlinear  relationship  between  stress 
and  displacement,  both  the  probability  density  function  and  the 
peak  distribution  of  displacement  response  process  are  signifi¬ 
cantly  different  from  those  of  the  linear  stress  response. 

The  number  of  stress  peaks  per  unit  time  and  the  up-urossing 
rates  increase  with  the  increase  of  input  sound  pressure  levels 
as  would  be  expected;  however,  the  nonlinear  stress  response  has 
a  mean  value  while  the  mean  value  for  the  nonlinear  displacement 
response  is  zero. 

The  spectral  densities  of  the  nonlinear  response  show  a  wid¬ 
ening  of  response  peaks  and  a  shift  towards  higher  frequencies  as 
the  input  levels  increase  and  the  nonlinearity  effects  become 
more  dominant. 

The  time  domain  analysis  presented  in  this  study  indicates 
that  for  anticipated  sound  pressure  levels  acting  on  present  and 
future  aircraft  structures,  the  various  simplified  linear  theo¬ 
ries  used  to  predict  stress  response  and  fatigue  life  would  not 
produce  realistic  structural  panel  configurations.  The  rather 
dramatic  differences  between  linear  and  nonlinear  predictions  is 
significant  and,  thus,  is  a  reminder  that  structures  exposed  to 
acoustic  noise  must  be  carefully  designed. 
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Figure  1  A  Rectangular  Panel  Exposed  to  Random  Pressure 
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Figure  2  Computer  Programs  Used  to  Perform  Time  Domain  Analysis  and  Predict  Fatigue  Life 


Spectml  density.  In  ■'Hz  Spectra!  density.  In  /Hz 


Figure  9  Spectral  Density  of  Displacement  Response  for  Input 
Sound  Pressure  Levels  of  80  dB,  100  dB,  and  120  dB 
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Figure  13  Probability  Density  and  Peak  Distribution  Histograms 
of  Displacement  w  for  150  dB  Input 
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Figure  14  Probability  Density  and  Peak  Distribution  Histograms 

of  Normal  Stress  a.,  for  120  dB  Input 
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Figure  15  Probability  Density  and  Peak  Distribution  Histograms 
of  Normal  Stress  a  for  150  dB  Input 
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igure  19  Displacement  Response  Time  History  (150  dB  Linear) 
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Figure  20  Displacement  Response  Time  History  (150  dB  Noni inear) 
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Figure  22  a  Stress  Response  Time  History  (150  dB  Nonlinear) 
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Probability  Density  Histogram  of  Displacement 
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Figure  28  Probability  Density  Histograin  of  o  (150  dB  Nonlitiear) 
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Figure  33  Up-Crossing  Rate  for  a  as  a  Function  of  Lay-Up  Angle 
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RMS  STRESSES  AS  A  FUNCTION  OF  ANGLE 
LAMINATE  I+THETA/-THETAlx2  SYMMETRIC 
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Figure  36  RMS  of  as  a  Function  of  Lay-Up  Angle.  C+0/-0]x2s 
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APPENDIX  A 


LISTING  OF  COMPUTER  PROGRAM  TDR 

C****  *♦** 

C****  TDR 

C****  *♦♦♦ 

c****  time  domain  response  analysis  of  a  simply  supported  plate  ***♦ 

C****  SUBJECTED  TO  UNIFORM  RANDOM  PRESSURE  AND  THERMAL  LOAD 

c#*** 

C****  BASED  ON  THE  WORK  OF  RIMAS  VAICAITIS  OF  COLUMBIA  UNIVERSITY 

AND  S.  T.  CHOI  (RESEARCH  ASSISTANT),  1988-89  ♦♦♦* 

C****  ♦♦♦♦ 

C****  MODIFIED  TO  INCLUDE  ORTHOTROPIC  MATERIALS  BY  ROCKY  ARNOLD, 

ANAMET  LABORATORIES,  INC.,  1989  ♦♦♦♦ 

C**** 

C**** 

C****  SIMPLY-SUPPORTED  SINGLE  PANEL  SUBJECTED  TO  UNIFORM  RANDOM 
C***»  PRESSURE  AND  THERMAL  LOADS. 

***» 

ct***  this  program  is  used  to  find  THE  DISPLACEMENT  AND  STRESS 
C***t  RESPONSE  TIME  HISTORIES  FOR  A  SIMPLY  SUPPORTED  PANEL  USING 
C***#  MODAL  ANALYSIS  IN  THE  TIME  DOMAIN  WITH  NMODEX  MODES  IN  THE 
C****  X-DIRECTlON  AND  NOMODFY  MODES  IN  THE  Y-DIRECTION. 

C****  **** 

c****  LOADING  IS  UNIFORM  AND  TEMPERATURE  DISTRIBUTION  IS  ASSUMED 

UNIFORM.  ♦♦♦♦ 

0*4:**  **** 

C****  PARAMETERS  ***♦ 

C****  NX  =  MAXIMUM  NO.  OF  MODES  IN  X-DIRECTION  **** 

C****  nY  -  MAXIMUM  NO.  OF  MODES  IN  Y-DIRECTION  **** 

c**»*  nXY  =  NX*  NY  **** 

0****  nTEQ  =  MAXIMUM  TOTAL  NUMBER  OP  EQUATIONS  (=:2*NXY)  ♦*** 

C****  NSTEPT  =  MAXIMUM  NO.  OP  TIME  STEPS  »*** 

C****  NMODEX  =  ACTUAL  NO.  OF  MODES  IN  X-DIRECTION  USED  IN  ANALYSIS  ♦*** 

0****  NMODEX  =  ACTUAL  NO.  OF  MODES  IN  Y-DIRECTION  USED  IN  ANALYSIS  **** 

C****  rmsD  =  ROOT  MEAN  SQUARE  OF  DISPLACEMENT  RESPONSE  **♦* 

C****  RMSX  =  ROOT  MEAN  SQUARE  OF  STRESS  RESPONSE  SIGMA(X)  **** 

C****  RMSY  =  ROOT  MEAN  SQUARE  OF  STRESS  RESPONSE  SIGMA(Y)  *♦** 

C****  rmSXY  =  ROOT  MEAN  SQUARE  OF  STRESS  RESPONSE  SIGMA(XY)  **** 

C****  rloAD  -  GAUSSIAN  RANDOM  PRESSURE  ♦*** 

C****  **** 

0:(c:|c4:************************************************************************* 

IMPLICIT  REAL*8  (A-H,0-Z) 

PARAMETER  (MXPARM=50,  NX=3,  NY=3,  NXYs:9,  NTEQ=18,  NSTEPT=:8192) 

DIMENSION  A ( NX , NY ) , PARAM ( MXPARM ) , Z ( NTEQ } , ZETAI J ( NX , NY ) , 

1  VIJMN(NX,NY,NX,NY),ZNTIJ(NX,NY),CO(NX,NY), 

2  WIJ(JX,NY),ZIJKLMNRS(NXY,NXY,NXY,NXY),QIJ(NX,NY), 

3  WIJ2(NX,NY) 

DIMENSION  YR(18) ,DYR(18) , Y1R(18) ,Y2R{ 18) ,Y3E{18) 
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REAL*4  TIME(IO) 

COMMON  /COMO/  XL,YL,V,ZETA11 
COMMON  /COMl/  QIJ,WIJ2,C0 
COMMON  /COM2/  XYL,I,J,K,L 

COMMON  /COM3/  NMODEX,NMODEY,NMODEXY,NDIMX,NDIMY,NDIMXY,ILIN 
COMMON  /COM4/  VIJMN.ZIJKLMNRS 
COMMON  /COM5/  IHEAT.ZNTIJ 

COMMON  /RNG/  NR,XR,YR,DYR,HH, JR, JMAX,HE,X0UT,IFBEQ,X1R,X2R,X3R, 

1  Y1R,Y2R,Y3R,T0L 

COMMON  /XFER/  ISTEP,DSTEP,DDT,RLOAD(8192 ) 

COMMON/  PROPS/  A11,A12,A22,A66,H,T1,T2,T3 
DATA  NDIMX,NDIMY,NDIMXY/3,3,9/ 

c****  read  in  input  data  ♦♦♦♦ 

Qift**1f*******t******^^t****i^***********t*t*t*t**t*******T^^:^t*t******t*******t* 

TIME(1)=SECNDS(0.) 

CALL  READ(X,Y,DSTEP,NSTEP,NEQ,PI,WIJ,PQ1,C1,C2,C3,C4, 

1  PX , PY , XYL , XL2 , YL2 , VIJMN , ZIJKLMNRS , IHEAT , ZNTIJ , 

2  STHERX,STHERY,SPL) 

TIME{2)=SECNDS{TIME(1))+TIME(1) 

NR=N£Q 

C******tt**iii4ii*t*t****t****t***ilii^t***t*i^itt*!¥t*t****t***t****t*t*****t**t****** 
C****  CALL  ROUTINE  TO  CALCULATE  TIME  DOMAIN  SIMULATION  OF  PRESSURE  ♦*** 

TIME(3)=SECNDS(TIME(2) )+TIME(2) 

CALL  SIMLOAD(SPL) 

TIME(4)i=SECNDS(TIME(3))+TIME(3) 

C****  INITIALIZE  S'JMMING  PARAMETERS 

G4i4!4t4!4:4e4e4:4!4t4(4i<t«4;:4e*4:4:4i:|c4:4c«4i«4:i)!*$)tc4c«i|t4c4;:i|!4!:|c:|t«4t4;*4c4e*«4:4ei|ii|e4;4c*i(c!4(«4E*4!4e*iit«*4t4c*4Ci(ci(i«4c4titE 

SUMD=0.0 

SUMD2=0.0 

SUMX=0.0 

SUMX2=0.0 

SUMY=0.0 

SUMY2=0.0 

SUMXY=0.0 

SUMXY2=0.0 

T=0.0 

XR=T 

C****  COMPUTE  DAMPING  TERM  ♦♦♦* 

G^t^^4:itc4e4:4i4:*4i9t:4i4‘4:«3|t4!3|(*)|E4:4t4Etitc4!4(44t4:4t4ii|t4e4t4ci4ti)c4c#4:>tE3|c4E4c4citc4c4(i|c*4EitE4t4‘ii(’(‘4:^4ti|t4‘^i|(4ti|c#4E^i|Ei|t^4: 

DO  5  I=l,NMODEX 
DO  5  J=l,NMODEY 

ZETAIJ(I,J)=ZETA11*(WIJ(1,1)/WIJ(I,J)) 

C0(I,J)=:2.*ZETAIJ(I,J)*W1J(I,J) 

WIJ2{I,J)=WIJ(I,J)*WIJ(I,J) 

5  CONTINUE 

C*  COMPUTE  PRESSURE  TERM 

Gi!  ):4t4:4‘4:)t:$$4e4c*4t4ciii9):4t4i)itt4:4i*4:4t^ite4!4i4i!4t:|c4!4t4!4c4!:^:^i):#:|ei4c!|(4c4e:(Ci(cs|c#:t::l(4E4t4(4e!t:4:#4E4c#4:4:4i)(^4e4c4i4c4c:t!4:4c*it! 
DO  20  I=l,NMODEX 
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DO  2v0  Jsl.NMODEY 
20  CONTINUE 

C*ilf**t*********>¥****if*t**********tt*-,*^*******t****tt*t********************** 
C****  SOLVE  PDEb  for  DISPLACEMENT 

Citit**t*********^f^*********^***********i^***********#*************************<* 

TIM£(5)sSECNDS(TIME(4)}+TIN£(4) 

DO  10  ISTEP=1,NSTEP 

TEND  =  FLOAT (I STEP )*DSTEP 
HH=DMI N 1 ( DSTEP , DDT ) *0 . 1 2  5 
XOUT=TEND 
CALL  RUNGE 
DO  30  I=l,NNODEX 
DO  30  J=l,NMODEy 
K=(I-l)*NMODEY+J 
A(I,J)=YR(K) 

30  CONTINUE 

C****  COMPUTE  STRESS  RESPONSE 

C4e4e4t*4:4i*4i4i:|!4!4ii(!4t*4:44:i|ti|‘4‘4!4titt4E4i4e4!t4:4!4:4t*4t4!i|ciiE4:iK4:i|eilc4c3|c44t4ti4c4c4t#4i«i|ti(c*$4e)|t!|:il(t*4e4ci|ti(e4c4tt4t4:ik* 

T=XR 

TD=DSTEP 

CALL  RESP0N(A,XL,YL,ISTEP,TD,C1,C2,C3,C4,PX,PY,XYL,IHEAT, 

1  STHERX , STHERY , SUMD , SUMD2 . SUMX , SUMX2 . SUMY , SUNY2 , 

2  SUMXY,SUMXY2) 

10  CONTINUE 

TIME(6)=SECNDS(TIME(5) )+TIME(5) 

C****  COMPUTE  MEAN,  MEAN  SQUARE  AND  RMS  VALUES 

(;;#4:4:4:4!4c#4!4:4:]|e4!#4t4!4c44c^#4!4‘#4:ii‘4(#*4i4*4:4i44!4!4c3ie4[4c4:)ic4c:tc4e4t4c)|e4(«4c#$4t4tiic4c>ic4!4c4t4t4:i|i44:i|c$t4c#<c4;<(4t4t 

SUMD=SUMD/FLOAT ( NSTEP ) 

SUMD2=SUMD2/FLOAT ( NSTEP) 

RMSD=DSQRT(SUND2) 

SUMX=SUMX/FLOAT ( NSTEP ) 

SUMX2=SUMX2/FLOAT ( NSTEP ) 

RNSX=DSQRT(SUNX2) 

SUMY=SUMY/FLOAT ( NSTEP) 

SUMY2=SUMY2/FLOAT ( NSTEP ) 

RMSY=DSQRT(SUMY2) 

SUMXY=SUMXY/FLOAT ( NSTEP ) 

SUMXY2=SUMXY2/FLOAT ( NSTEP ) 

RMSXY=DSQRT ( SUNXY2 ) 

Ctt*t  WRITE  OUT  SIMPLE  STATISTICS  OF  TIME  DOMAIN  RESPONSE  *♦** 

WRITE (6, 1000)  SUMD,  SUMD2,  RMSD 
WRITE(6,1003)  SUMX,  SUMX2,  RMSX 
WRITE(6,1004)  SUMY,  SUMY2,  RMSY 
c  WRITE(6,1005)  SUMXY,  SUMXY2,  RMSXY 

1000  FORMAT(’  Displ.  (in):  Mean  =  ’ ,E11.4, ’  M.S.=  ’,E11.4, 

+  ’  RMS  =  ’,E11.4) 

1003  FORMAT(’  sigaaX  (psi):  Mean  =  ’,E11.4,*  M.S.=  ’,E11.4, 

+  ’  RMS  =  *,E11.4) 


1004  FORMAT(’  sigBaY  (psi):  Mean  =  M.S.=  ’,E11.4, 

+  ’  RMS  =  *,E11.4) 

clOOS  FORMATS  tauXY  (psi):  Mean  =  ’ ,E11. 4, ’  M.S.=  *,E11.4, 

c  +  *  RMS  =  *,E11.4) 

C** 

PRINT  1100 

1100  FORMAT(/i*  Output  files  FOR008:  Response  histories',//) 

TIME ( 7 ) =SECNDS ( TIME ( 6 ) ) +TIME ( 6 ) 

C****  PRINT  OUT  TIME  SUMMARY  *♦♦♦ 

WRITE(6,9999) 

9999  FORMAT (///,X, 'SUMMARY  OF  TIME  EXPENDITURES*,/) 

WRITE(6,9998)  TIME(2)“TIME{1) 

9998  FORMAT{/,X,’TIME  IN  READ  SUBROUTINE=  *,F8.1,’  SECONDS') 

WRITE(6,9997)  TIME{3)-TIME(2) 

9997  FORMAT(/,X,*TIME  BETWEEN  READ  AND  SIMLOAD  SUBROUTINES^  *,F8.1, 

1  *  SECONDS') 

WRITE(6,9996)  TIME(4)-TIME(3) 

9996  FORMAT (/,X, 'TIME  IN  SIMLOAD  SUBROUTINES  ',F8.1,*  SECONDS') 

WRITE(6,9995)  TIME(5)-TIME(4) 

9995  FOBMAT(/,X,‘TIME  BETWEEN  SIMLOAD  AND  RUNGE  SUBROUTINES=  *,F8.1, 

1  '  SECONDS') 

WRITE(6,9994)  TIME(6)-TIME(5) 

9994  FORMAT{/,X,'TIME  IN  RUNGE  SUBROUTINES  *,F8.1,*  SECONDS’) 

WRITF(6,9993)  TIME(7 )-TIME(6) 

9993  FORMAT!/, X, ’TIME  BETWEEN  RUNGE  AND  FROGRAM  END=  *,F8.1, 

1  ’  SECONDS’) 


**** 

C**t*  I  J  K  L  ***^'* 

*4t4t4: 

C**1litt*******t*t***t***^tt***ilf***t*t**t**********tt***^****ii****tt*tt******** 
SUBROUTI NE  I JKLMNRS ( VIJMN , Z 1 JKLMNSS , ZNTIJ , P2 , P3 , P4X , P4Y , 

1  XL2,YL2,IHEAT) 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  /COM2/  XYL,I,J,K,L 

COMMON  /COM3/  HMODEX , NMODEY . NMODEXY , NDIMX , NDIMY , NDIMXY , I LI N 

COMMON  /PROPS/  A11,A12,A22,A66,H,T1,T2,T3 

DIMENSION  VIJMN(NDIMX, NDIMY, NDIMX, NDIMY) ,ZNTIJ(NDIMX, NDIMY) , 

1  Z I JKLMNRS ( NDIMXY , NDIMXY , NDIMXY , NDIMXY ) 

c****  COMPUTE  LINEAR  (HOMOGENEOUS)  COMPONENT  OF  AIRY  STRESS  FUNCTION  *♦♦♦ 

DO  10  Isl,NMODEX 
CI2XL2=I«I/XL2 
DO  10  Jsl, NMODEY 
CJ2YL2=J*JAL2 
DO  10  M-l,NMODEX 
CM2XL2:=M«M/XL2 
DO  10  Nsl, NMODEY 
CN2YL2=N*NAL2 

VIJMN( I , J , M , N ) ^P2* ( Cl 2XL2* (T1*CM2XL2+T2*CN2YL2 ) 

1  +CJ2YL2* (T2*CM2XL2+T3*CN2YL2 ) ) 

10  CONTINUE 

IFdLIN  .EQ.  0)  GO  TO  20 

C****  COMPUTE  NONLINEAR  (P/ITICULAR)  COMPONENT  OF  AIRY  STRES  FUNCTION  »**♦ 
Qilf****tt***t***************  .:.***rlll*****ift********************tt**************** 
DO  15  I=l,NNODEX 
DO  16  J=l, NMODEY 
IJ=(I-l)*NMODEY+J 
DO  15  K=l,NMODEX 
DO  15  L=l, NMODEY 
KL=(K-l)*NMODEY+L 
DO  15  M=l,NMODEX 
DO  15  N=l, NMODEY 
MN=(M-l)*NMODEY+N 
DO  15  IR=l,NMODEX 
MPRsM-JlR 
MMR=M-IR 

DO  15  1S=1, NMODEY 
IRS=(IR-l)*NMODEY+IS 
NPS=N+IS 
NMS=N-IS 
NR=N*IR 
MS-M*IS 

ZIJKLMNRS(IJ,KL,MN,IRS)  =  P3  « 

1  (MS*(NR-MS)*(FIJ(MPR,NPS)+FIJBAR(MMR,NMS))+ 

2  MS*(NR+MS)*(FIJ{MMR,NPS)+FIJ(MPR,NMS))) 

Id  CONTINUE 
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c«««*  COMPUTE  THERMAL  COliPONENT  OF  AIRY  STRESS  FUNCTION  **** 

20  IF  (IHEAT  .NE.  1)  THEN 
ELSE 

DO  30  I=l,NMODEX 
CI2XL2=I*I/XL2 
DO  30  J=l,NMODEY 

ZNTIJ  ( I ,  J )  =P4X*CI  2XL2+P4Y*  J*  J  AL2 
30  CONTINUE 
END  IF 
RETURN 
END 
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C****  **** 
C****  S  I  M  L  0  A  D  **♦♦ 
C^*!***  **♦* 


C 

c 

c 

c 

c 

c 


N  —  MO.  OF  INTERVALS  IN  THE  SPECTRUM 

N  SHOULD  BE  AN  INTEGER  POWER  OF  TWO 
NPT  —  NO.  OF  POINTS  FOR  THE  TIME  SERIES 
NPT  SHOULD  BE  INTEGER  POWER  OF  TWO. 
ISEED  ~  RANDOM  NUMBER  SEED 


NPT>N 


SUBROUTINE  SIMLOAD(SPL) 

IMPLICIT  BEAL*8  (A-H,0-Z) 

COMMON  /XFER/  ISTEP,DSTEP,DT,Y(8192) 

DIMENSION  X(8192) ,SP(3S00) ,W(3500) ,BAND(8192) 

COMPLEX  X,ZIMAG 
LOGICAL  INVERSE 

DATA  FMAX,INVEBSE/500.,.TRUE./ 

DATA  N.NPT  /S12,  4096/ 

C4:«4cik**)|i*»4!*4!«**********»««««««4t««««4t3|t»««»**«!(!4c4t4t4t***4c*4>****iii**4‘***iic4c**«***** 


C****  INITIALIZE  VARIABLES  *♦♦♦ 

SPP=8.41*10**(-18.+SPL/10. ) 

PI  »  3.141592654 

PI2  s  PI  *  2.0 

NPl  =  N  1 

Z1MA6  s  CNPLX(0. 0,1.0) 

SPPW=SPP/PI2 

WUaFMAX*PI2 

DW  -  WU  /  FLOAT(N) 

DO  119  Isl,NPl 
SP(I)=SPPW 
W(I)=:(I-1)*DW 
119  CONTINUE 

AREA:kSPP*FNAX 
SQ2DW  s  DSQRT(2.0«DW) 

TTOTAL=PI2/DW 
Dt-ttOTAL/FLOAT ( NPT ) 

C****  SET  X(1)=0.  IN  ORDER  TO  OBTAIN  NEW  MEAN  ZERO  TIME  SERIES  *♦** 

X(l)  =  CMPLX(0. 0,0.0) 

DO  50  I  =  N+1,NPT 
X(I)  =  CMPL.^(0. 0,0.0) 

50  CONTINUE 


C****  GENERATE  RANDOM  PHASE  ANGLES  UNIFORMLY  DISTRIBUED  BETWEEN  **** 

C****  ZERO  AND  2.*PI  **** 

C:(tt#4c4t4t44c)):!|(44‘#4c#4‘i|t##4e4t4‘i|ti|t3|t4E)|t4c$4t4t4t)|‘4i**4c4c##4(4t4t4i4:4t4t4i4>4E4:)it4t4t!ii4‘4ti|(#4t4t4i*4:4‘4:4t*#4:4c4!|t4‘4c4‘ 

ISEED=12357 
DO  51  Isl,N 

51  RANlCnsRANdSEED) 

DO  60  1=2, N+1 


105 


PHI  =  RAND(I-l)  *  PI 2 
PI  =  SQ2DK  ♦  DSQBT(SP(I)) 

X(I)  =  PI  *  CDEXP(-ZIMAG*PHI) 

60  CONTINUE 

c*«**  PERFORM  FORWARD  TRANSFORM  **** 

CALL  FFT  (X,NPTa) 

get  real  PART  ♦♦♦*  ^ 

DO  70  1=1, NPT 
Y(I)  =  REAL(X(I)) 

70  CONTINUE 
RETURN 
END 
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C«*«***»«*«*«*««*««*«*«4i««»»«»«*«*»«**4t**«****»«»»*****»*»*4i«»»4t»«»»*»«*«¥»*« 
C**t*  «*»« 

C«*«*  F  I  J  **** 

C****  ’ll*** 

FUNCTION  FIJdQ.IH) 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON  /CON2/  XYL,I,J,K,L 

COMMON  /PROPS/  A11,A12,A22,A66,H.T1,T2,T3 

IPKal+K 

IMK-I-K 

JPLsJ+L 

JMLsJ-L 

KHsK«IH 

LGsL*IG 

FIJs2 .  «KH«LG*  ( BETA  (I  PK ,  I G ) -t-BETA  ( I  NX ,  I G ) )  *  ( BETA  ( JPL ,  I H  )  i 

1  BETA(JNL,IH))-(KH««2-l-LG**2)*(6AMMA(IPK,IG)- 

2  GAMMAdNK, IG)  )»(GAMMA(mdH)-GAMMA(  JML,  IH) ) 
DEN0M=A22*IG**4+(A66+2.*A12)*XYL**2*<IH»IG)**2 

1  +A11*XYL**4*IH**4 

FIJ=FIJ/(H*DENOM) 

RETURN 

END 


**** 

C««4:*  F  I  J  B  A  R  **** 

C*4ct«  *«t« 

C*:^*t*****t**il^********t****t*\‘*t*******t****t**tif********t******************* 

FUNCTION  FIjMR(IG,IH) 

IMPLICIT  RMAL^e  (A-H,0-Z) 

COMMON  /COMV  XYL,I,J,K,L 
COMMON  /PROPS/  A11,A12,A22,A66,H,T1,T2,T3 
IF  (IG  .EQ.  IH  .i\ND.  IG  .EQ.  0)  THEN 
FIJBARsO.O 
ELSE 
IPK=I+K 
IMK=I-K. 

JPL=J+L 

JMLsJ-L 

KH=K*IH 

LG:sL«IG 

FIJBAR=2.*KH*LG*(BETA(IPK,IG)+BETA(IMK,IG))*(BETA(JPL,IH)+ 

1  BETA(JML,IH))-(KH**2+LG*»2)*<GAMMA(IPK,IG)- 

2  GAHMA(IMK,IG))«(GANMA(JPL,IH)-GANNA(JML,IH)) 
DEK0M=A22*IG**4+(A66+2.*A12)*XYL**2*(IH*IQ)**2+A11*XYL**4*IH**4 
FIJBARsFIJBAR/ (H*DENOM ) 

END  IF 
RETURN 
END 


C****  **** 

c****  BETA  **** 

C**t*  **** 

C»**«»*«**«*««*»»»»«««»»««*«*»*»«**««i|t*****««««*»«*»»**««**««*«*#**«»**«»«*«« 

FUNCTION  BETAdP.IQ) 

IMPLICIT  BEAL«8  (A-H,0-Z) 

IF  (IP  .EQ.  IQ  .AND.  IP  .NE.  0)  THEN 
BETA  s  1.0 

ELSE  IF  (IP  .EQ.  -IQ  .AND.  IP  .NE.  0)  THEN 
BETA  =  -1.0 
ELSE 

BETA  =  0.0 
END  IF 
RETURN 


€«**♦  **** 

C****  GAMMA  ♦♦♦♦ 

C*»*»  **** 

FUNCTION  GAMMA(IP,IQ) 

IMPLICIT  BEAL*8  (A-H,0-Z) 

IF  (IP**2+IQ**2  .EQ.  0)  THEN 
GAMMA  a  2.0 

ELSE  IF  (lABS(IP)  .EQ.  lABS(IQ)  .AND.  IP  .NE.  0)  THEN 
GAMMA  =  1.0 
ELSE 

GAMMA  -  0.0 
END  IF 
RETURN 
END 


C««*»«««»«««»«»*«*«***4t«**)|i4t*«»**««*«*«**««****«*«««*«*«**«*««*»«*»««*«*«t4c’l‘« 
C****  *♦** 

D  1  P  F  E  Q  ****‘ 

C****  **** 

C*««««*«***»««»*»»«»*»»««»«»*««»*««*»*««»»***»*******««*«4i«*«**)K*4:««««*«**««« 
SUBROUTINE  DIFFEQ 
IMPLICIT  BEAL«8  (A-H,0-Z) 

COMMON  /RNG/  N,X.Z,ZPRIME,HH,JR,JNAX,MR.XOUT,IFB£Q, 

1  X1,X2,X3,Y1,Y2,Y3,T0L 

COMMON  /XFES/  I3TEP,DSTEP,DT,ELOAD{8192) 

DIMENSION  Z(18),ZPRIME(18).Y1(18),Y2(18),Y3(18) 

DIMENSION  QIJ(3,3),WIJ2(3,£hOO(3,3),ZNTIJ(3,3), 

1  VIJMN(a,3,3,3).r]LJKLMNSS(9,9,9,9) 

amm  /com/  qij,  wij2,  cu 

COMMON  /COM3/  NMODEX,NMODEY,NMODEXY,NDIMX,NDIMY,NDIMXY, ILIN 
COMMON  /COM4/  VIJMN,  ZIJKLMNBS 
COMMON  ZOOMS/  IHEAT.  ZNTIJ 
DATA  ICNT/0/ 

C«»««  INTERPOUTE  TO  DETERMINE  LOAD  TERM  **** 

NEQ=N 

IF(ICNT.GT.O)GO  TO  2 

ICNT=1 

I PLUS* 1 

PRMsO.O 

PRPaO.O 

TlsO.O 

TOa-DT 

SLOPEaO.O 

PRaO.O 

2  iF(X.GT.Tl)  GO  TO  1 

IF(X.LT.TO)  GO  TO  3 

PSaPRM+(X-TO)*SLOPE 
GO  TO  20 

3  IPLUS=“1 

1  ICNTaICNT+IPLUS 

PRMaBLOAD(ICNT-l) 

PRPaRLOAD(ICNT) 

SLOPE* (PRP-PRM)/DT 
T0aT0+DT*IPLUS 
TlaTl+DT*IPLUS 
I PLUS* 1 
GO  TO  2 
20  CONTINUE 

C4t«««  SPECIFY  DIFFERENTIAL  EQUATIONS 

DO  5  Kal.NMODEXY 
ZPRIM£(K)sZ(K-t-NMODEXY) 

5  CONTINUE 

DO  10  lal.NMODEX 
DO  10  Jal,NMODEY 


K=(I-1)»NM0DEY+J 

KK^NNODEXY-I-K 

CALL  VZIJ(VIJ,ZIJ.ZtNEQ,I,J) 

ZPe!NE(KK)  =  PB*QIJ(I.J)-C0(I,J)*Z(KK)-WIJ2(I,J)«Z(K) 
1  -VIJ*Z(K)-ZIJ+ZNTIJ(I,J)*Z(K) 

10  CONTINUE 
RETURN 
END 
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C**tt  **** 

c***t  V  Z  I  J  **** 

C4t*»«  *«»* 

SUBROUTINE  VZIJ(VIJ,ZIJ,Z,NEQ,I , J) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  Z(NEQ),  VIJNN(3,3,3.3) , 

4-  ZIJKLMNRS(9,9,9,9) 

COMMON  /COM3/  NMODEX,NMODEY,NMOD£XY,NDIMX,NDIMY,ND!MXY,ILIN 
COMMON  /COM4/  VIJNN,ZIJKLNNRS 

C** 

VIJsO.O 

ZIJsO.O 

IJ=(I-l)#NMODEY+J 

C** 

DO  10  M=l,NMOD£X 
DO  10  Nsl.NMODEY 
MN=(M-l)*NMODEY+N 
IFdLIN  .EQ.  0)  GO  TO  10 
VIJ=VIJ+Z  (MN )  *3  ( MN ) ‘VUMN  ( I ,  J ,  M ,  N ) 

ZKLRSsO.O 
DO  20  K=l,NMODEX 
DO  20  L=l,NMODEY 
KL=(K-l)*NMODEY+L 
ZRS^O.O 

DO  30  IRsl.NMODEX 
DO  30  ISsl.NMODEY 
IRS=(IB-l)*NMODEY+IS 
ZBS=ZRS-I-Z  ( I RS )  «Z  I JKLNNRS  ( IJ ,  KL ,  MN ,  I RS ) 

30  CONTINUE 

ZKLRS=ZKLRS-I-Z  ( KL }  *ZRS 
20  CONTINUE 

ZIJ£ZIJ-l-Z(MN)«ZKLRS 
10  CONTINUE 

RETURN 
END 
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C«**«  «««* 

C****  E  E  S  P  0  N  »*♦* 

4:*** 

C***************ift***t****i)(t*t***iti*!¥*********t****1ft**t*t*******t******t**tt* 
SUBROUTINE  BESPON ( A , XL , YL , I STEP , DT , Cl , C2 , C3 , C4 , PX , PY , XYL , I HEAT , 

1  STHERX , STHERY , SUMD , SUND2 , SUNX , SUMX2 , SUMY , SUMY2 , 

2  SUMXY,SUMXY2) 

IMPLICIT  REALMS  (A-H,0-Z) 

DIMENSION  A(NDIMX,NDIMY) 

COMMON  /COM3/  NMODEX,NMODEY,NMODEXY,NDIMX,NDIMY,NDIMXY,ILIN 
COMMON  /PROPS/  A11,A12,A22,A66,H,T1,T2,T3 
C4c4:*4c4t4!**4i*4t4c4:**i|:*««4i«***«4i4E«4c4i4t4:4:4E4c««««4c*4t4e4e4t*»i(t4c*4Ei|i»»iK!(t«4i*«4!4:4e4Ei|c4e4::^4c4ciK«!tc4c«$ 

c****  initialize 

C4e4(4c**)<t*4i**4t4:4t4ti|t4!4!4t4i*«*«4:4i4!4!4!*4e«4:4t4i*««4i**4c4c4t4c4c*4e4:*«4:4c4t***4e*4(4t4:«4(«*4::^*4:4t*i|c4i«4:* 

DISPL=0.0 

SIGMAXsO.O 

SIGMAY^O.O 

TAUXY=0.0 

C«4:«)|!*«»**4i«4E«**4:««4!4!4i4i4!«4i4c4!4c*4!«4:««4!4:4t«4t««i|c»4e4c*«*«»4t*4e«*4t*4t4c4c***i|E*4:4c4c«4t4c4iiiCiic4:A 
C****  COMPUTE  STRESSES  BY  SUMMING  LINEAR  AND  NONLINEAR  TERMS 
(;;4e4t44:4!4!i|!4i4t4!4:4E««4i4:4:4t4:4:«**4i4i*4t*«4<4E4c4!«4t*4‘4t4t4c«4t4:«4t«i|eiit*«***4c4t»«»««*4c4:4(«*t4(«4i4e4(*4:»» 
DO  10  M=l,NMODEX 
SINMX=DSIN(M*PX) 

COSMX=DCOS(M*PX) 

XLM=M/XL 

DO  10  N=l,NMOD£Y 
YLNsNAL 
SINNY=DSIN(N*PY) 

COSNYsDCOS(N*PY) 

DISPLsDISPL+A(M,N)*SINMX*SINNY 

C*«**  CALCULATE  STRESSES  FROM  NONLINEAR  PART  OF  AIRY  STRESS  FUNCTION  **** 

CALL  SUMRS ( RSXl , RSX2 , RSYl , RSY2 , RSXYl , RSXY2 , M , N , 

1  A,PX,PY,XL,YL,XYL) 

CC1=:C2*A(M,N) 

IF(ILIN  .EQ.  0)  CC1=0. 

C*4i«4!4c4t«*4:*«4:*4i«4c4t4<4:4t4c4t«**«4:4!««4!«»****4t4c*4:4c*«*«**«4c*4(4i4:4:**4c4c4i4«4ci|t4E«4ii(t*4i4t)|t«*4c* 
C«*«*  ADD  LINEAR  AND  NONLINEAR  CONTRIBUTIONS 

C4:4t4i4!«*4cs|i4E4c**4:4(4i4!4!**4c4c4!4!#4t«4c«4E4c4t4i4:4t4:4c4t4c*4c*«*4(**»**!t:*4(4i*4c»*4(»«4c4(«**4c4(4c4E**it:«4(«« 
CC0=A (M , N ) * (C1*SI NMX*SI NN . +CC1 ) 

SIGMAX=SIGMAX+CC0* (T1*XLM**2+T2*YLN**2 ) -C3*A(M , N) ♦ (RSX1+RSX2 ) 
SIGMAY=SIGMAY+CCO* (T2#XLM**2+T3*YLN**2 ) -C3*A(M, N)* (RSY1+ESY2 ) 
TAUXY=TAUXY+ ( C4/A66 ) ♦ A ( M , N ) *COSMX*COSNY-C3*A ( M , N ) ♦ ( RSXYl +BSXY2 ) 

10  CONTINUE 

C4i**4(4c«4e4t4!*4e**4!4E4i»»4E4e«)|!***4c)^4i4E4!4t»4i*«4(4e«4(i|e*4t»««i|t4i4c«4t«)lt«4c!|E4e4t4t««*««*«««4t*4(4!«««*i(i 
C«4:*»  add  in  THERMAL  COMPONENT  **** 

0*4: 4!»4i*4:«*4(«4c4E4(*«4e««4c4c4c4i*4i  4(4(4(41 4:  4e*4e  4c  4c  4c«4t4E4(*4c*4c4ti|c4!4t4Ei|:4!4:4i*4t*««)|c«»4i4i4:4t4:*4e««4!4e)|t4(  4(4: 
IF  (IHEAT  .NE.  1}  THEN 
ELSE 

SIGMAX-SIGNAX'STHERX 
SI GMAY=SI GMAY-STHERY 


114 


END  IF 

C»***  PRINT  OUT  TIME  DOMAIN  BESPONSE  **** 

T=DT*ISTEP 

NBITE(6,100)  T,  DISPL,  SIGNAX,  SIGMAY 
100  FOBMAT(2X,F8.6,3(2X,E14.6)) 

WRITEd.lOOO)  T, DISPL, SIGNAX, SIGMAY, TAUXY 
1000  FORMAT (5E10. 3) 

C«**«  SUM  DISPLACEHENT/STBESSES  AND  SQUARES  *»*» 

SUNDsSUND’I'DISPL 

SUND2sSUMD2-l-DISPL«DISPL 

SUNXaSUMX+SIGMAX 

SUMX2sSUNX2-«-SIGNAX«SIGMAX 

SUNY=SUMY-f  SIGMAY 

SUNY2sSUMY2-»-SIGMAY«SIGNAY 

SUNXYsSUMXY4-TAUXY 

SUMXY2=SUMXY2+TAUXY»>TAUXY 

RETURN 

END 
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*««« 

C****  S  U  N  R  S  **** 

C***«  **** 

SUBROUTINE  SUMRS ( RSX 1 , RSX2 , RSYl . RSY2 . RSX Y1 , RSX Y2 , N , N , 

1  A,PX,PY,XL,YL,XYL) 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON  /COM3/  NMODEX,NNODEY,NMODEXY,ND1MX,NDIMY,NOIHXY,ILIN 
COMMON  /PROPS/  A11,A12,A22,A66.H,T1,T2,T3 
DIMENSION  A(NDIMX,NDIMY) 

C4:*i|E4:4t4(4tM#4i)ki|t4:4c*^3|c4t)t:4i4t4!i|>3|t*44t4:4ii|t4i4i4!«i|i*4t4t*4:4:4i>lt4t4i|t4t4ci|t4()|E4Ei|:4t4t:(ti|t4!!|i4ti|c!(c4c4ttts(iiii)|t4t«4:4t*4s 

C****  INITIALIZE 
S1=A22*H 

S2=(2.*A12iA66)«H 

S3=:A11*H 

RSXl-0.0 

RSX2-0.0 

RSY1=0.0 

RSY2s0.0 

RSXYlsO.O 

RSXY2=0.0 

IF(ILIN  .EQ.  0)  RETURN 

COMPUTE  NONLINEAR  STRESSES  (PROM  PARTICULAR  SOLUTION)  **** 

DO  10  IRsl.NMODEX 
MPRsM+IH 
MMR=M-IR 
MPR2aMPR*MPR 
MMR2=MMR»MMR 
XLMPR=MPR/XL 
XLMMR=NNR/XL 
XLMPR2-XLMPR««2 
XLMMR2sXLMMR*«2 
CMPRXsDCOS(MFR«PX) 

CMMRXsDCOS(MMR*PX) 

SMPRX=DSIN(MPR*PX) 

3NMRX=:DSIN(MMR«PX} 

DO  10  IS=l,NMODEY 
NPS=N-fIS 
NMS=N-IS 
NPS2=NPS*NPS 
NMS2=NMS*NMS 
YLNPS=NPS/YL 

ylnms=nmsal 

YLNPS2=YLNPS*»2 

YLNMS2=YLNMS**2 

CNPSY=DCOS(NPS*PY) 

CNMSY=DCOS(NMS*PY) 

SNPSY=DSIN(NPS*Py) 

SNMSY=DSIN(NMS*PY) 
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XyUPS=XYL*NPS 

XYNMSbXYL«NNS 

XYNPS2*XYNPS**2 

XyNM32aXYNM3**2 

C#«*  TERMS  FOR  WHICH  "KR  .NE.  ML” 

IE  (N*IE  .EQ.  M«IS)  THEM 
ELSE 

PPt=CMPRX*CNPPV/(31*MPR2*«2+S2*MPR2*XYNPS2+S3*XYNPS2»*2) 
PP.*.-CNNRX«CNMIiY/(Sl«MMR2«*2+S2»KMR2«XYNNS2+S3»XYNNS2**2) 
PP5-SIIPRX«SNPSY/ ( Sl*MPR2««2-l>S2«MPR2*XYMPS2+S3«XYNPS2««2 ) 
PP6sSNMRX«SNMSY/(Sl«MHR2*«24S2«MMR2«XYNMS2+S3*XYNMS2»*2) 
Q1»A(IB,IS)«N«IS«(H*IR-N*1S) 

RSXlsRSXl'fQl«(  YLN[PS2«PP14-YLRMS2«PP2  ) 

RSY1»R3Y1-Hai«  ( XLNPR2«PP1-I-XLNKR2*PP2  ) 

RSXYlsR3XYl-V<21«  (  XLMPR«YLNP3«PP5^XLXMR«YLNMS«PP6 ) 

END  IE 

C***#  TEEMS  FOR  WHICH  ”KR  .EQ.  ML”  ***■* 

PP3sCNPBX*CNMSY/(Sl«HPR2«»2+S2«MPR2«XYNMS2-l-S3*XYNHS2**  { ) 
PP4aCMMRX*CNPSY/ ( S1*MMB2**2+S2»MMR2*XYNPS2+S3*XYNPS2**  >. ) 
PP7aSMPRX*3»:<4&Y/(Sl*MPa2**2+S2CMPB2#XYNMS2+S3*XYNMS2**2) 
PP8sSMMEX«SNPSY/  ( Sl«MME2««24-S2*NMB2*XYNPS2-i'S3»XYNPS2**2 ) 
Q2sA(lR,IS)»N*IS«(N*IR^H«IS) 

RSX2sBSX2+Q2* ( YLNMS2*PP3+YLNPS2*PP4 ) 

BSY2sBSY2-t-Q2«  ( XLMPB2*PP3-t-XLNME2*PP4 ) 

RSXY2sRSXY2+Q2* (XLMPB*YLNMS*PP7+XLMMR*YLNPS#PP8 ) 

10  CONTINUE 
RETURN 
END 
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C**** 

C****  R  U  N  G  E  **** 

SUBROUTINE  RUNGE 
IMPLICIT  BEAL«8  (A-H,0-Z) 

COMMON  /BNG/N,X,Y,DY,HH,J,JMAX,M,  X0UT,IFBEQ,X1,X2,X3,Y1,Y2,Y3,T0L 
DIMENSION  Y(18),DY(18),Y1(18),Y2(18),Y3(18) 

J=1 

JMAX^l 

IFREQ=3 

M=1 

CALL  EUNKUT 

RETURN 

END 
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C«***«*««*»«»*»«»»«»«»«*««***«««***»***»«***«*«**»»»**»***«*««««*«**«««*t.«»*« 

c****  «**« 

C****  R  U  N  K  U  T  ♦♦** 

C**** 

SUBROUTINE  RUNKU'i' 

IMPLICIT  R£AL*8  (A-H,0-Z) 

COMMON/RNG/N,  X,Y,DY,HH, J, JMAX.M,  X0UT,IFREQ,X1,X2,X3,Y1,Y2,Y3,T0L 
DIMENSION  Y(18),DY(18),Y1(18),Y2(18),Y3(18) 

INDE9  -  0 

CALL  ADJSTP 

IF{J-JMAX)  10,10,50 
10  INDE9  »  INDE9  +  1 
CALL  INTPOL 
IF(J<JMAX}  20,20,50 
20  CALL  STEP 
XI  =  X2 

X2  =  X3 

X3  =  X 

DO  30  I  =  1,  N 
Y1(I)  =  Y2(I) 

Y2(I)  =  Y3(I) 

30  Y3(I)  s  Y(I) 

IF  (INDE9  >  IFREQ)  10,40,40 
40  INDE9  s  0 

CALL  ADJSTP 

IF(J-JNAX)  10,10,50 
50  RETURN 
END 
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c***» 

c**** 

ADJSTP 

0**** 

C********************t*********^***********t******i*t*t********ftfilittt********* 
SUBROUTINE  ADJSTP 
IMPLICIT  REAL*8  (A“H,0-Z) 

COMMON/RNG/N,X,Y,DY,HH,J,JMAX,M,XOUT,IFREQ,Xl,X2,X3,Yl,Y2,Y3,TOL 
DIMENSION  Y(18) ,DY(18) ,Y1(18) ,Y2(18) ,Y3(18) 

KSL=0 

HFACT  =  1.0  D+31 

HFACTl  =  l.OD+30 

GO  TO  (30,10),  M 
10  HI  =  HH 

HH  =  2.0  ♦  HH 

X  =  XI 

DO  20  I  -  1,  N 

20  Y(I)  =  Y1(I) 

GO  TO  100 
30  KSL=1 

40  HI  HH 

XXX  =  X 

DO  50  I  =  1,  N 
50  Y1(I)  =  Y(I) 

XI  =  X 

CALL  INTPOL 
IF(J-JMAX)  60,60,250 

60  CALL  STEP 
DO  70  I  =  1,  N 
70  Y2(I)  =  Y(I) 

X2  *  X 

CALL  INTPOL 
IF(J-JMAX)  80,80,250 
80  CALL  STEP 

DO  90  I  =  1,  N  . 

Y3(I)  =  Y(I) 

80  Y(I)  =  Y1(I) 

X3  =  X 

X  =  XXX 

HH  -  2.0  HH 

100  CAL^  STEP 

DO  150  I  =  1,  N 

DELY  =  DABS  (  Y(I)-Y3(I ) )/30.0 
IF(DELY  -DABS  (Y2(I))*TOL  )120,110,110 
110  IF(  DABS  (Y2(I))-TOL)  120,130,130 
120  HFIRST  l.OD-t-30 
GO  TO  140 

130  HFIRST=  (DABS  (Y2(I))*  TOL/DELY  )  **0.2 
140  CONTINUE 

150  HFACTsDMINl  (HFACT,  HFIRST  ) 

IF  (HFACTl  -  HFACT)  160,160,170 
160  HH  =  2.0  *  HI 

GO  TO  (40,230),  H 
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170  HH  =  HI  ♦  HFACT 

GO  TO  (180,230),  M 
180  IF(KSL)  220,220,190 
190  KSL=0 

IF(DABS  (HH)-DABS  (HI))  200,220,220 
200  DO  210  I  -  1,  N 


210  Y(I) 

=  Y1(I) 

X 

-  XXX 

GO  TO  40 

220  KSL-0 

M 

=  2 

230  DO  240  I 

=  1,  N 

240  Y(I) 

-  Y3(I) 

250  RETURN 

END 

C******»***********#i|t**#***********<!***t*j((*****)|c:*!****i|.*«j|tKt)(t*#****j),**4t**»***** 

C**t* 

C****  STEP 

C****  4i4!4c4t 

C’^****t*****t*******illi*ttttt^*t*t****t*t**t***tt**tt**************t*********** 

SUBROUTINE  STEP 
IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  /RNG/N,X,Y,DY,HH,J,JMAX,M,XOUT,IFREQ,Xl,X2,X3,Yl,y2,Y3,TOL 
DIMENSION  Y(18),DY(18),Y1(18),Y2(18),Y3(18) 

DIMENSION  Y0(18),P1{18) 

DO  10  I  =  1,  N 
10  Y0(I)  =  Y(I) 

XO  =X 

CALL  DIFFEQ 
DO  20  I  =  1,  N 
P1(I)  =  DY(I)  ♦  HH 

20  Y(I)  =  Y0(I)  +  P1(I)*0.5 

X  =  XO  +  HH*0.5 

CALL  DIFFEQ 
DO  30  I  =  1,  N 
P1(I)=  P1(I)+2.0*HH*DY(I) 

30  Y(I)  =  Y0(I)  +  0.5*HH*DY(I) 

CALL  DIFFEQ 

DO  40  I  =  1,  N 

P1(I)=  P1(I)+2.0*HH*DY(I) 

40  Y{I)  =  Y0(I)  +  HH*DY(I) 

X  =  XO  +  HH 

CALL  DIFFEQ 

DO  50  I  =  1,  N 

50  Y(I)=YO(I)  +  (P1(I)+HH*DY(I))^0.1666667 
RETURN 
END 
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c**** 

c**** 

INTPOL 

tt** 

c**** 

SUBROUTINE  INTPOL 
IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  /RNG/N,X,Y,DY,HH,J,JMAX,M,X0UT,IFREQ,X1,X2,X3,Y1,Y2,Y3,T0L 
DIMENSION  Y(18),DY(18),Y1(18),Y2(18),Y3(18) 

IF{DABS  (XOUT  -  X)-DABS  (HH))  10,10,20 
10  HH=XOUT-X 
CALL  STEP 
J  =  J  +  1 
20  RETURN 
END 


C****  *♦*♦ 

C*»»»  F  F  T 

Ct**if*****t***1ti*t$*****t******t****t**t*******>¥***t*ft:iit*******t:tt**t**illiM*tt*** 
SUBROUTINE  FFT{X,N,K) 

IMPLICIT  INTEGER  (A-Z) 

REAL«4  GAIN,PI2,ANG,BE,IM 
COMPLEX  X(N),XTEMP,T,U(16),V,W 
LOGICAL  NEW 

DATA  PI 2. GAIN, NO, KO/6. 283185307,1. 0,0,0/ 

NEW-NO.NE.N 
IF(.NOT.NEW)GO  TO  2 
L2N=0 
N0=1 

1  L2N=L2N+1 

N0=N0+N0 

IF(NO.LT.N)  GO  TO  1 
GAIN=1.0/N 
ANG=PI2*GAIN 
RE=COS(ANG) 

IM=SIN(ANG) 

2  IF(.NOT.NEW.AND.K*KO.GE.l)  GO  TO  4 

U{1)=CMPLX(RE,-SIGN(IM,FL0AT(K))) 

DO  3  Is2,L2N 

3  U(I)=U(I-1)*U(I-1) 

KO=K 

4  SBY2sN 

DO  7  STAGE=1,L2N 
V=U( STAGE) 

W=(l. 0,0,0) 

S=SBY2 

SBY2=S/2 

DO  6  L=1,SBY2 

DO  5  1=1, N,S 

P=I+L-1 

Q=P+SBY2 

T=X(P)+X(Q) 

X(Q)=(X(P)-X(Q))*W 

5  X(P)=T 

6  W=W*V 

7  CONTINUE 

DO  S  1=1, N 

INDEX=1-1 

JNDEX=0 

DO  8  J=1,L2N 

JNDEX=JNDEX+JNDEX 

ITEMP=INDEX/2 

I F  ( I TENP-^  I  TEMP .  NE .  1 NDEX )  JNDEX= JNDEX-i- 1 
INDEXsITEMP 

8  CONTINUE 

J=JNDEX+1 
IF(J.LT.I)GO  TO  9 
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XTENPsX(J) 

X(J)=X(I) 

X(I)=XTENP 

9  CONTINUE 

IF(K.6T.0)RETURN 
DO  10  I-1,N 

10  X(I)*X(I)*GAIN 

RETURN 

END 


PBOGRAM  PDF 


C***«  ***« 

c**»«  THIS  PSOQRAN  IS  USED  TO  CALCULATE  THE  PROBABILITY  DENSITY  «*«« 

C««»«  FUNCTION,  PEAK  DISTRIBUTION,  AND  UP-CROSSING  RATE  OF  A  RANDOM  **«* 
C**»*  PROCESS.  **** 

C****  **** 

c**«*  the  CALCULATION  IS  DONE  FROM  -4*SD  TO  4*SD  WITH  NDIV  INTERVALS 
C**««  IN  EACH  STANDARD  DEVIATION  (SD).  THE  TOTAL  NUMBER  OP  INTERVALS  «««* 
C**t*  IS  8*NDIV  WITH  THE  FIRST  INTERVAL  BEING  -INFINITY  TO 
C*«*»  -{4,-(4-l/NDIV)SD)  AND  THE  LAST  INTERVAL  BEING  (4-l/NDIV)SD)  TO 

C****  INFINITY. 

C«**«  «*«* 

c«*»«  the  input  data  file  is  F0RQ01.DAT  WHICH  CONTAINS  THE  RESPONSE 
C«»««  TINE  HISTORY  PRODUCED  BY  PROGRAM  TDR.  **** 

if*** 


c**************************************************************************** 
COMMON/TITLES/  TITLE(20) ,SUBTIT(20) 

DIMENSION  F(13000),  XDENS(200),  XNPK(200),  XCROSS(200) 

CALL  REAODATA( F , NPT , SD , DSD , NHALF, NHALFl , 

1  NDIV,NDIVT,DX,TTOTAL,SHIFT) 

CALL  DENSITY ( F , XDENS , NPT , SD , DSD , NHALF , NHALFl , 

1  NDIVT,DX,SHIFT) 

CALL  PEAKDI S ( F , XNPK , NPT , SD , DSD , NHALF , NHALFl , 

1  NDIVT,DX,TTOTAL,SHIFT) 

CALL  UPCROSSR ( F , XCROSS , NPT , SD , DSD , NHALF , NHALFl , 

1  NDIVT,DX,TTOTAL, SHIFT) 

C««*»***4«*«**«**»«««««»»««*«**«*«««»«4t***««i»*««»«$»««*»*»»«»****«**«*«*««*«4i 

**** 

c****  INFORM  USER  OF  OUTPUT  FILES  AND  CONTENTS 

C****  **** 

Qt*************************************************************************** 
TYPE  ’ 

TYPE  Output  files:  FOR007  —  Probability  Hensity  and’ 

TYPE  theoretical  Gaussian’ 

TYPE  ♦,*  Output  files:  FOR008  —  Peak  distribution’ 

TYPE  4,*  Output  files:  FOR009  —  Upcrossing  rate* 

TYPE  *•’  Output  files:  FOROlO  —  Theoretical  Gaussian’ 

TYPE  ♦,’  Output  files:  FOROll  —  Theoretical  Rayleigh’ 

STOP 

END 


C****  SUBROUTINE  READDATA  ♦♦♦♦ 
c****  **** 
C****  READS  DATA,  COMPUTES  THE  MEAN,  RNS,  STANDARD  DEVIATION,  **** 
C****  COEFFICIENTS  OF  SKEWNESS,  KURTOSIS  OF  THE  PROCESS  **** 
C**** 


SUBROUTINE  READDATA ( F , NPT , SD , DSD , NHALF , NHALFl , NDI V , NDI VT , DX , 
1  TTOTAL,SHIFT) 

DIMENSION  F(l) 

COMMON/TITLES/  TITLE(20) ,SUBTIT(20) 
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CHAR4CTER«8G  DUMP 

C«»««  READ  IN  DATA  FROM  TDR  (FOR001.DAT)  AND  SCREEN  (UNIT  5)  **** 

REAO( 1,10000)  TITLE 

10000  FORNAT(X,20A4) 

BEAD( 1,10000)  SUBTIT 
READd, 10001)  NPT,DT 

10001  FORMAT(I5,E10.3) 

READ  (5,111)  ICOL,NmV,IMEAN 
111  FORMAT(3I5) 

TTOTALsDT*(NPT-l) 

DX=:1./NDIV  !  DX  =  nondisensional  increment. 

NDIVT=NDIV*8  !  Total  #  of  divisions  from  -4SD  to  +4SD 

NHALF=4*NDIV 

NHALFlsNHALF+1 

IF  (ICOL  .EQ.  2)  THEN 

DO  81  1=1, NPT 

81  READ  (1,99)  ZJUNK,  F(I) 

99  FORMAT (5E10. 3) 

ELSE  IF  (ICOL  .EQ.  3)  THEN 

DO  82  1=1, NPT 

82  READ  (1,99)  ZJUNK,  ZJUNK,  F(I) 

ELSE  IF  (ICOL  .EQ.  4)  THEN 

DO  83  1=1, NPT 

83  READ  (1,99)  ZJUNK,  ZJUNK,  ZJUNK,  F(I) 

ELSE 

DO  84  1=1, NPT 

84  READ  (1,99)  ZJUNK,  ZJUNK,  ZJUNK,  ZJUNK,  F(I) 

END  IF 

find  mean  and  RMS  FOR  PROCESS  **** 

C**t**t*^******if^*****************t***lft****************if*t$**********!¥***’lt** 
SUN=0.0 
SUM2=0.0 
DC  10  1=1, NPT 
SUM=SUM+F(I) 

SUM2=SUM2+F(I)**2 
10  CONTINUE 

XMEAN=SUM/NPT 

SUM2=SUM2/NPT 

RMS=SQRT(SUM2) 

C**t*  FIND  STANDARD  DEVIATION  (SD),  COEFFICIENTS  OF  SKEWNESS,  KURTOSIS  *♦** 

SUMV=0.0 

SUMS=0.0 

SUMK=0.0 

DO  20  1=1, NPT 

DIFF=F(I)-XMEAN 

DIFF2=DIFF*DIFF 

SUMV=SUMV+DIFF2 

SUM'3=SUMS+DI  FF*DI  FF2 


SUMKsSUMK+DI FF2«DI FF2 
20  CONTINUE 

SD=SQBT(SUMV/NPr) 

COSKEW=SUMS/SD**3/NPT 
COKURT=SUMK/SD=t'*4/NPT-3 . 

C*»«*»***»*««**««*««***«4i»»**«»**»»*«t»«»»***»4c«3i!««»««««*»«««*««»»»*«*««*«*** 

C****  WRITE  OUT  CALCULATIONS 

WRITE(6,1111) 

nil  FORMAT(1H1,//,10X,’  PDF 

1  lOX,’  FOR  AFSC/ASD/PMRNA 

2  lOX,’  WPAFB,  OHIO 

3  lOX,*  BY  ANAMET  LABORATORIES,  INC.*,/, 

4  lOX,*  HAYWARD,  CALIFORNIA  ’,/, 

5  lOX, ’CONTRACT  NO.  F33615-89-C-3210’ ,/) 

WRITE(6,1234)  (TITLE(KK) ,KK=1,20) 

1234  FORMAT (/,X,20A4) 

WRITE(6,1234)  (SUBTIT(KK) ,KK=1,20) 

WRITE  (6,30)  XMEAN,  RMS,  SD,  COSKEW,  COKURT 
30  FORMAT(10X,’  MEAN  =  ’ ,E13.5,/, 

+  lOX,*  RMS  =  ’,E13.5,/, 

+  lOX,’  STANDARD  DEVIATION  =  *,E13.5,/, 

+  lOX, ’COEFFICIENT  OF  SKEWNESS  =  ’,E13.5,/, 

+  lOX, ’COEFFICIENT  OF  KURTOSIS  =  ’,E13.5,/) 

C****  DIAGNOSTIC  MESSAGE  *♦♦♦ 

SHIFT=XMEAN/SD 

IF  (SHIFT  .GT.  0.05)  THEN 

WRITE  (6,100)  XMEAN, SD, SHI FT 

100  FORMAT  (/,’  MEAN  VALUE  OF  THE  PROCESS  IS  NOT  ZERO!!!!  ’,/, 

1  ’  MEAN  s  ’,E13.5,’  SD=’,E13.5,’  RATIO  s  * ,F7.4,/) 

ELSE 
END  IF 

C««««  RECOMPUTE  F(I)  WHEN  MEAN  IS  EXCLUDED  **** 

DO  98  1=1, NPT 
F(I)=F(I)-XMEAN 
98  CONTINUE 

IF  (IMEAN  .EQ.  1)  THEN 

SHIFT=0.0 

ELSE 

END  IF 

DSD=SD/NDIV  !  Increment  in  actual  value 

C****  WRITE  RESULTS  ON  OUTPUT  FILES 
WRITE  (7,310) 

310  FORMAT(’  PROBABILITY  DENSITY’,/) 

WRITE  (7,313)  XMEAN, SD, RMS, NPT,DT,NDIVT 
WRITE  (7,314)  IMEAN 
WRITE  (8,311) 


311  FOHNAT(*  PEAK  DISTBIBUTION*) 

WRITE  (8,313)  XN£AN,SD,RMS,NPT,DT,NDIVT 
WRITE  (8,314)  INEAN 
WRITE  (9,312) 

312  FORMAT (’  UP-CROSSING  RATE*,/) 

WRITE  (9,313)  XM£AN,SO,RNS,NPT,DT,NDIVT 
WRITE  (9,314)  IMEAN 
WRITE  (10,315)  SD 
WRITE  (11,316)  SO 

313  FORMAT  (*  MEAN  =  *,E13.5,*  SD  -  *,E13.5,*  RMS  s  *,E13.5,/, 

1  ’  NO.  OF  PTS  =’,I6,*  DT  =  ’,E13.5,’  TOTAL  DIV.  *  *,I3) 

314  FORMAT(/,2X,*  IMEAN  a:  *,I2,*  NOTE:  IF  IMEAN-1  THEN  MEAN  EXCLUDED’) 

315  FOSMAT(*  THEORETICAL  GAUSSIAN  WITH  MEAN  -  0.  AND  *, 

1  ’  SD  s:  ’,E13. 5, /////) 

316  FORMAT( ’  THEORETICAL  RAYLEIGH  WITH  * , 

1  ’  SD  =  ’,E13. 5, /////) 

RETURN 

END 

Ct***  SUBROUTINE  DENSITY  **** 

€♦**♦  THIS  SUBROUTINE  CALCULATES  THE  PROBABILITY  DENSITY,  THEORETICAL  **** 

C****  GAUSSIAN  DESNTIY  WITH  ZERO  MEAN  AND  STANDARD  DEVIATION  OF  THE  **** 

C****  RANDOM  PROCESS  *♦** 

SUBROUTI N£  DENSITY ( F , XDENS , NPT , SD , DSD , NHALF , NHALFl , NDI VT , DX , SHI FT ) 
COMMON/TITLES/  TITLE(20) ,SUBTIT(20) 

DIMENSION  F(l),  XDENS(l) 

PI-3.1415926 

SQ2PI-SGRT(2.*PI) 

DO  10  1=1, NDI VT 
XDENS(I)=:0. 

10  CONTINUE 

DO  20  1=1, NPT 
TEMP=F(I)/DSD 

ITEMP=IINT(TEMP)  !  Locate  the  data  belongs  to  which  interval 

IF  (TEMP  .GE.  0.)  THEN 

IF  OTEMP  .GT.  NHALF)  ITEMP=NHALF*1 

XDENS ( NHALF1+ I TEMP ) =XDENS ( NHALFl+I TEMP ) + 1 . 

ELSE 

IF  (I TEMP  .LT.  -NHALF)  ITEMP=-NHALF+1 
XDENS  ( NHALF-f  I  TEMP )  =XDENS  ( NHALF+ITEMP )  -•■  1 . 

END  IF 

20  CONTINUE 

c****  WRITE  RESULTS 

WRITE(7,1234)  (TITLE(KK) ,KK=1,20) 

1234  FORMAT(X,20A4) 

WRITE(7,1234)  (SUBTIT(KK) ,KK=1,20) 

WRITE  (7,210) 

210  FORMATS  MAGNITUDE/SD  PROBA.  DENS.  NO. OP  OCCUR.  GAUSSIAN’, 

1  /,’  (NORMALIZED)’) 

WRITE  (10,211) 


211  FORNAT('  HAGNITUDE/SD  PBOBA.  DENSITY*, 

1  /,’  (GAUSSIAN)’) 

AREA=0.0 

X=-4.0 

DO  40  I=1,NDIVT 
XDX=X+DX 

XDENS1=XDENS( I )/NPT/DX 
ABEA=:AREA-t-XDENSl 

C****  CALCULATE  THE  THEORETICAL  GAUSSIAN  WITH  ZERO  MEAN  **** 

GAUSS=EXP(-X*X/2. )/SQ2PI 
GAUDX=EXP(-XDX*XDX/2. )/SQ2PI 
YsX+SHIFT 
YDX=Y+DX 

IF  (I  .EQ.  1)  WRITE  (7,220)  Y, ZERO, ZERO, GAUSS 
WRITE  (7,220)  Y,  XOENSl,XDENSO  )  fOAUSS 
WRITE  (7,220)  YDX,XDENS1,XDENS(I ) ,GAUDX 
220  FORMAT  (2(3X,F11.5),3X,F8.1,3X,F11.5) 

IF  (I  .£Q.  NDIVT)  WRITE  (7,220)  YDX, ZERO, ZERO, GAUDX 

IF  (I  .EQ.  1)  WRITE  (10,220)  Y,  GAUSS 

WRITE  (10,220)  Y,  GAUSS 

WRITE  (10,220)  YDX,  GAUDX 

IF  (I  .EQ.  NDIVT)  WRITE  (10,220)  YDX,  GAUDX 

XsX+DX 

40  CONTINUE 

AREAsABEA«DX 

WRITE  (6,700)  SD,  AREA 

700  FORMAT  (*  SD  =  ’,E12.5,’  AREA  OF  DENSITY  CURVE  *  *,F8.4,/) 

RETURN 

END 


C**tt  **** 

Ct***  SUBROUTINE  PEAKDIS 

C*«»«  THIS  SUBROUTINE  CALCUUTES  THE  PEAK  DISTRIBUTION  AND  **** 

C««**  the  THEORETICAL  RAYLEIGH  DISTRIBUTION  *♦♦♦ 


SUBROUTI NE  PEAKDI 3 ( F , XNPK , NFf , SD . DSD , NHALF , NHALFl , NDI VT , DX , 

1  TTOTAL, SHIFT) 

COMMON/TITLES/  TITLE(20)  ,SUBTItUo) 

DIMENSION  F(1),XNPK(1) 

DO  15  1=1, NDIVT  ^ 

XNPK(I)=0. 

15  CONTINUE 

DO  20  1=1, NPr  t 

IF  (I  .EQ.  1  .OR.  I  .EQ.  NPT)  GO  TO  L 

DF1=F(I)-F(I-1) 

DF2=r(I+l)-F(I) 

DSIGN=DF1*0F2 

IF  (DSIGN  .GT.  0.)  GO  TO  20 
IF  (DF2  ,GT.  PFl)  GO  TO  20 
TEMP=F(I)/DSD 
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20 


30 


1234 


210 


211 


212 


ITEMPsnNT(TEMP) 

IF  (TEMP  .GT.  0.)  THEN 

IF  (ITEMP  .GT.  NHALF)  ITEMP^NHALF-l 

XNPK(NHALFl-l>ITEMP)sXNPK(NHALFl4>ITEMP)-l-l. 

ELSE 

IF  (ITEMP  .LT.  -NHALF)  ITEMP=-NHALF+1 
XNPK(  NHALF+ITENP)  =XNPK(  NHALF+ITENP)-!-! . 
END  IF 
CONTINUE 


XNPEAK=0. 


DO  30  I=1,NDIVT 
XNPEAK=XNP£AK-i-XNPK(I ) 

CONTINUE  !  XNPEAK  IS  THE  TOTAL  NO.  OF  PEAKS 

PEAKTsXNPEAK/TTOTAL  !  NO.  OF  PEAKS  PER  UNIT  TIME 

WHITE  (6.210)  XNPEAK.  PEAXT 
WRITE(8,1234)  (T1TL£(KK) .KK=1.20) 

FORMAT(X,20A4) 

WRITE(8.1234)  (SUBTIT(KK) ,KKsl,20) 

WRITE  (8.210)  XNPEAK.  PEAKT 

FORMAK’  TOTAL  NO.  OF  PEAKS  =  ‘.FlO.l, 

1  *  NO.  OF  PEAKS  PER  SEC.  =  *,F10.1) 

WRITE  (8.211) 

PORMAT(*  MAQNITUDE/SD  PEAK  DI STB.  NO.  OF  PEAKS  RAYLEIGH*. 

1  /,*  (NORMALIZED)*) 


WRITE  (11.212) 

FORNAT(*  MAGNITUDE/SD  PROBA.  DENSITY* 
1  /.’  (RAYLEIGH)*) 


Xs-4.0 


DO  40  lal.NDIVT 
XDXs:X-t-DX 

XNPKlaXNPK( I )/XNP£AK/DX 
IF  (XDX  .LT.  0.0)  THEN 
RAY=0.0 


RAYDX=0.0 


ELSE 

RAY=X*EXP(-X*X/2. ) 
RAYDX=XDX*EXP(-XDX*XDX/2. ) 
END  IF 


Y=X+SHIFT 


YDX=Y+DX 

IF  (I  .EQ.  1)  WRITE  (8,220)  Y, ZERO, ZERO, RAY 
WRITE  (8,220)  Y,  XNPKl,XNPK(I ) ,RAY 
WRITE  (8,220)  YDX,XNPK1,XNPK(I),RAYDX 
220  FORMAT  (2(3X,F11.S) ,3X,F8.1,3X,F11.5) 

IF  (I  .EQ.  NDIVT)  WRITE  (8,220)  YDX, ZERO, ZERO, RAYDX 

IF  (I  .EQ.  1)  WRITE  (11,220)  X,  RAY 

WRITE  (11,220)  X,  BAY 

WRITE  (11,220)  XDX,  RAYDX 

IF  (I  .EQ.  NDIVT)  WRITE  (11,220)  XDX,  RAYDX 

X=X+DX 

40  CONTINUE 
RETURN 
END 
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4til|t4c4c 

SUBKOUTINE  UPCBOSSB 

«««* 

THIS  SUBROUTINE  CALCULATES  THE  UPCE03SINQ  >?ATE  AT  DIFFERENT  **«* 

THRESHOLD  LEVELS.  *♦♦♦ 

*»»* 

eUBROUTI NE  UPCROSSR ( F , XCROSS , MPT , SD , DSD , KHALF , NHALFl , NDI VT , DX , 

TTOTAL.SKIFT)  ^ 

COMMON/TITLES/  TITLE(20) ,SUBTIT(20) 

DIMENFTON  F(l),  XCBOSS(l)  v 

mTE{ 6,1234)  (TITLE(KK),KK=1,20) 

FOBMAT{X,20A4) 

?VRITE(9,1234)  (SUBTIT(KK),KK-1,20) 
mTE  (9,210) 

FORMAT ( / , 3X , ’ THRESHOLD  LE VEL/3D  UPCROSSI NG  RATE  ( #/SEC . ) ' ) 

DO  10  I=1,NDIVT 
ICBOSS(I)sO. 

CJONTINUE 
DO  20  1=1, NPT 

IF  (F(I)  .GT.  F(I+1))  GO  TO  20  !  Only  up-crosaing  is  counted. 

TEMP1=F(I)/D3D 

TEMP2=F(I+1)/DSD 

ITEMPl=nNT(TEMPl) 

ITEMP2=IINT(TEMP2) 

IF  (TEMPI  .GT.  0.)  THEN 

DO  30  K=ITEMP1,ITEMP2-1 

XCROSS  ( NHALFl4l-t-K)  =XCROSS(  NHALFl-t-l+K)  ‘t-1 . 0 

CONTINUE 

ELSE  . 

END  IF 

IF  (TEMP2  .LT.  0.)  THEN 

DO  40  K=ITEMP1,ITEMP2-1 

XCROSS (NHALF+l+K ) =XCBOSS ( NHALF+l+K) +1 . 0 

CONTINUE 

ELSE 

END  IF 

IF  (TEMP1*TEMP2  .LT.  0.)  THEN 
DO  50  K=I TEMPI, 0 

XCBOSS(NHALF+1+K)=XCROSS(NHALF+1+K)+1.0 

CONTINUE 

DO  60  K=0,ITEMP2-1  t 

KCROSS(NHALF1+1+K)=XCROSS(NHALF1+1+K)+1,0 

CONTINUE 

ELSE  ^ 

END  IF 

CONTINUE 

^=-4.04DX  !  the  lowest  level  for  crossing  is  -(4-DX)SD 

30  70  1=1, NDI VT 

if~X+SHIPT 

SCROS3 ( I ) =XCROSS ( I ) /TTOTAL 
JiaiTE  (9,JuIl)  Y,  XCROSSd) 
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200  FORMAT  (10X,F11.5,10X,F9.3) 
X=X+DX 

70  CONTINUE 
RETURN 
END 


C****  **** 

Ct***  This  progrsM  calculates  the  dasage  fros  a  given  peak  **** 

C****  distribution.  Input  the  standard  deviation  of  the  stress 

C*«a4r  process,  SD  (in  ksi)  and  the  total  #  of  peaks  per  second, 

C;*«4e4!  TPEAK,  and  the  aaterial  fatigue  constants  ZLANDA  and  B. 

Xhe  peak  distribution  histogran  is  found  on  FOB008  created  **t* 
cane**  froa  the  prograa  PDF  run  earlier  (10  divisions  per  s.d.  in 

the  histograii).  Output  is  E[Nt]*tau  and  tine  to  failure  in  **** 
C*»*»  seconds. .  *♦** 

C***t 


WRITE(6,1010} 

1010  FORMAT (1H1,//,10X,’  DAMAGE  *,/, 

1  lOX,’  FOB  AFSC/ASD/PMENA  *,/, 

2  lOX,’  WPAFB.OHIO  * ,/, 

3  lOX,*  BY  AKAMET  LABORATORIES,  INC.*,/, 

4  lOX,*  HAYWARD,  CALIFORNIA  *,/, 

5  lOX, ’CONTRACT  NO.  F33615-8S-C-3210’ ,//) 
fiEAD(S,100)  SD, TPEAK, ZLAMDA.B 

100  FORMAT (4F10.0) 

WRITE(6,1009)  SD, TPEAK, ZLAMDA,B 

1009  FORMAT (lOX,*  STANDARD  DEVIATION  OF  STRESS  PROCESS^  *,F10.1,/, 

1  lOX,*  TOTAL  NUMBER  OF  PEAKS  PER  SECOND>  ’,F10.1,/, 

2  lOX,*  FATIGUE  PARAMETER,  LAMDA^  *,F10.2,/, 

3  lOX,*  FATIGUE  PARAMETER,  B=  *,E10. 3,/) 

DXaO.ltSD 

SUMsO.O 
READ  (8,1100) 

1100  FORMAT (/////////) 

C** 

DO  10  1=1,80 
READ  (8,1200)  X,Y 
1200  FORMAT(2G,/) 

IF  (Y  .EQ.  0.0)  GO  TO  10 

X=(X+0.05)*SD 

PR0B=Y/SD 

SLAMDA=ABS ( X ) ««ZLAMDA 
SUN=SUM+PBOB«SLAMDA 
10  CONTINUE 
SUM=SUM«DX 
EMT=B/SUM 
TFAIL=B/SUM/TPEAK 
WRITE  (6,1300)  SUM,  ENT,  TFAIL 
1300  FORMAT(10X,’  INTEGRATION  =  ’,E14.6,/, 

1  lOX,*  S[Mt]*tau  =  ’,E14.6,/, 

2  lOX,’  TIME  TO  FAILURE=  ’,Ei4.6,’  SEC’,/) 

STOP 

END 


<rU.S.  Ckivamnwht  Printing  0((ic«:1»92— 648-127/S2638 
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DEPARTMENT  OF  THE  AIR  FORCE 

MU  FORCE  RESEARCH  LABORATORY 
WRiQHT-PArrensoN  air  force  base  ohio  45433 


MEMORANDUM  FOR: 


Defense  Tedmical  Infonnation  Center/OMI 
8725  John  J.  Kingman  Rd,  Suite  0944 
Ft  Belvoir,  VA  220606218 


FROM:  Det  1 AFRL/WST 

Bldg  640  Rm  60 
2331 12th  Street 

Wright-Patterson  AFB  OH  45433-7950 

SUBJECT;  Notice  of  Changes  in  Technical  Rqjortfs)  .  ( ^ 
Please  change  subject  teport(s)  as  follows: 


i(ji^pd-r/?-9o  -30?/) 
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